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THE INCREMENTAL PHASE PLANE 
FOR NONLINEAR SAMPLED DATA SYSTEMS 


John A. Aseltine and Richard A. Nesbit 
Space Technology Laboratories, Inc. 
University of California 
Los Angeles, California 


Abstract 


An analysis technique for nonlinear sampled 
data systems is developed, using the incremental 
phase plane. This method is analogous to the 
phase plane method for continuous systems. A 
sampled data system with saturation is analyzed 
to demonstrate the use of the incremental phase 
plane in system analysis. Path tangent curvesare 
introduced which allow the graphical solution of 
systems with more general types of nonlinearities. 
A difference equation analog of Van der Pol's 
equation is solved using the path tangent curves. 


Introduction 


Difference equations bear the same relation 
to sampled data systems as differential equations 
do to continuous ones. This suggests that there 
might be a counterpart to the phase plane, but 
based on differences, which could be applied to 
nonlinear sampled-data systems. This paper 
describes such a technique by introducing the 
incremental phase plane or A-plane. Inthe A- 
plane we plot the first difference (as opposed to 
the first derivative) against the system variable. 
It can be applied to sampled systems in a way 
similar in some respects to the application of 
phase plane to continuous systems. 


Before beginning the development of the 
method it is wellto callattentionto other methods 
proposed for nonlinear sampled systems. The 
conventional phase plane was first used by Alter 
and Helstrom! to analyze oscillatory systems. 
Kalman“ uses a "normal mode" phase plane in 
which system variables are plotted after a coor- 
dinate transformation. More recently Mullin3 has 
combined conventional phase plane with switching 
lines indicating sampling times. Schmitt? treats 
the nonlinearity as an equivalent gain in his 
studies. The first three of these methods are 
based on the conventional phase plane. That is, 
they make use of derivatives, rather than differ- 
ences, in their plots. 


In this paper we shall begin by showing the 
genesis of the difference equation describing the 
sampled system. Next, the A-plane will be intro- 
duced and A-plane trajectory equations defined 
and obtained. As an example, a system with 
saturation will be analyzed--an example of 
interest since the other graphical methods? 3 
have had difficulty treating it. Finally path tan- 

ent curves will be introduced and used to study 
a Van der Pol type sampled system. 


The A-plane method is applicable to second- 
order systems containing a variety of nonlinear- 
ities. In addition to saturation discussed in this 
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paper, the method has been applied to systems 
with bang-bang and other types of control for a 
number of plant transfer functions.2 We now 
proceed to develop the steps leading to the 
A-plane. 


Genesis of the Difference Equation 


Although the difference equationis the natural 
expression of behavior of a sampled system, we 
arrive at it in a round about way. This is due in 
part to the way in which sampled-data analysis 
evolved. We will here summarize the steps 
briefly. 


Figure 1 shows a linear system described 
by transfer function G(s). At its input and output 
are samplers which produce impulse functions 
periodically, each having an area equal to the 
value of the sampled variable at that instant. 


FIGURE |. SAMPLED LINEAR SYSTEM 


The samplers are shown by the switches in the 

figure. For convenience we shall normalize the 
sampling period to unity. An input is applied at 
time zero. The sampled input, then, is given by 


x sia 
en(t) = E  e,(t) 6(t-n) (1) 
n=o 
The X transform of this is 


co 
EZ e,(n) eo (2) 
n=o 


e;(s) = 


If we now replace 2 by z we have the A - 
transformed input 


e;(n) gin a [es | (3) 


Note that in (3) we have also defined the % 
transform. 


(ee) 
E,(z) = 5 
TX=O0 


See reference 6 for details of sampled systems, 
and reference 7, Chap. 16 for % transforms. 


Next, the {,-transformed output is given by 


E,(s) G(s) E*(s) (4) 


the derivation 


Without going into the details qf 
we have for the 


which can be found elsewhere, 
%-transform of the output 


E, (2) G(z) E,(z) (5) 


where 


G(z) = (6) 


{x mol 


* 

in which h is the sampled impulse response. We 
now have means of expressing the transforms of 
sampled input and output of a linear system. 


Let us turn now to a question concerning the 
4%, transform defined by (3). If we were to begin 
a table of operation-pairs, we might start with 


e(n) 


e(n+1) 


where the question arises regarding the operation 
in the z-domain corresponding to shifting in n. 
From (3) we write 


[6 @) 
Bos) e(nats ae 7) 
n=oO0 


fore) 
Zle(nt1)] = EX e(nti)z™ 


n=0 


Expanding the sum on the right we have 


Pee ier Lely co ct e(2)e Ke. 


iT) 


z[E(z) - e(0)] (8) 


and we can fill in the table 


E(z) 


e(n+1) zE(z) - ze(0) 

We see that multiplication by z corresponds to 
shifting in n in the same way that multiplication 
of a Laplace transform by s corresponds to 
differentiation. 


Returning now to (5) we write 


E (z) a EIT ho a 

eg) > Ge) - 444+ = (9) 
i bee ai been oon ated b 
a r-1 fo) 


where we have assumed that G(z) is a rational 


fraction. Multiplying (9) out, we have 
r r-1 
(b,. Zz + Lalo Z +. + b,) E, (2) 
2 q ak 
(Ayz rs Eze eet: A.) E,(z) (10) 
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and, finally, using the shifting property (8) and 
assuming initial conditions are included in the 
input, the difference equation results: 


be (mtr) =f: bly e,(ntr-t) cere tars, bo e(n) 
(11) 
= RY) + ee ae ae? to eet A, e; (n) 
Thus, the correspondence 
e(n+1) Zt (Z))eeee (0) ee O (12) 
and the one related to the first difference 
Ae(n) (z-1) E(z) , e(0) = 0 (£3) 
where 
Ae(n) = e(nti) - e(n) = (E-1t)e (14) 


become the keys to the difference-equation des- 
cription of a sampled system. 


As an example, consider the double integrator 
shown in Figure 2. 


52 


eke Le 


ej 


FIGURE 2.DOUBLE INTEGRATOR 


The z block diagram is shown in Figure 3 
with the corresponding z transfer function. 


Balle 


FIGURE 3..z BLOCK DIAGRAM 
OF DOUBLE INTEGRATOR 


Now we write 


Oe AE een: (15) 
iL (z-1) 

or 
(21)° Ej(z) = 2, (2) (16) 


and, using (12) and (13), the system difference 
equation results: 


AS e,(n) = 


e(n+1) (17) 


The A-Plane 


The incremental phase plane, or A-plane, is 
shown in Figure 4. Because ofthe definition ofthe 
A operator (14) the ordinate at any point is equal 
to the difference in abscissa to the next sample 
point. 


FIGURE 4. A—PLANE 


Each point of a solution specifies xn and 
Axn. The value of xn}{ is obtained from xp+4 = 
x, + Ax,- The value of Ax,44 can be determined 
by either of two methods. One method uses tra- 


jectories, the other uses path tangent curves. 


A trajectory is defined to be a curve which 
passes through all points of a solution that has 
any point on the curve. 


It should be pointed out that the values of 
sampled system variables between sampling in- 
stants are not defined along the trajectory, which 
only connects values at the sampling instants. 
When an ambiguity arises, as between points A 
and B inthe figure, the furthest point reached 
along the trajectory is the correct one. 


Trajectories 


Let us write the second-order difference 
equation 


Ane =K (18) 
as two first-order equations 
Ay = K 
(19) 
Ax = y 
The equation of the trajectory for the system. 
described by these equations is 
y* - K(2x+y) = C (20) 


as can be verified by differencing (20). [ A method 
of direct solution of (19) in the form of (20)--that 
is, without first knowing the solution to the paired 
equation2--is not known to the authors.] The 
trajectory (20) is given by a family of parabolas 
offset from the x axis, one of which is shown in 
Figure 5. 


AX 


FIGURE 5, TRAJECTORY FOR A®x =| 


Another case of interest is the equation for a 
purely oscillatory system with frequency w, 


2 
A x + 2(1-cos w,) Ax + 2(1-cos wo) x =) — (2%) 


which we write 


Ay + K(y+x) = 0, K © 2(t-cosu) (22) 
The trajectory equation is 
2 
oa Vesa ts eC ee (23) 


a family of rotated elipses as shown in Figure 6 
for K=1. Equation (23) can be verified as a 
solution of (22) by differencing. 


Ax 


FIGURE 6. TRAJECTORIES FOR A*+AX+X=0 


These trajectories will be used in the example 
treated in the next section. 


A System with Saturation 


We shall illustrate the application of the A- 
plane method by considering the system with 
saturation shown in Figure 7. 
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FIGURE 7. SYSTEM WITH SATURATION 


From (17) we write (assuming the input zero) 


Rela) = Safer) (24) 
But x, the output of the limiter has a value de- 


pending on c, so 


1 c(n+1) <-1 
Nectatusedece(n et) Jc(nt1)| <4 (25) 
-1 c(n+1) >1 


In other words, if the subsequent value of c lies 
between +1, the equation is 


A clalenotiesy) 266 (26) 


or from (14) 


A@c(n) + Ac(n) + c(n) = 0 (27) 


and the trajectory is given by Figure 6. 


On the other hand, if the subsequent value of 
c lies outside the saturation zone, the trajectory 
is given by Figure 5 or its inverted image, since 
then 

NS eS (28) 

To construct a phase portrait of the system be- 
havior, we start with an initial point c(o), Ac(o), 
and use the appropriate trajectory to find sub- 
sequent points. At each point we determine 
whether the following one lies within a different 
region as defined by (25). If a new region is 
entered, the trajectory is changed accordingly. 


Several phase portraits for this system are 
shown in Figure 8. For the initial conditions 
chosen, limit cycles of periods 9, 10, and 14 
result. The dased parts of the portraits corres- 
pond to trajectories chosen from Figure 6. The 
other parts come from Figure 5 or its inverted 
image. 
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FIGURE 8.PHASE PORTRAITS FOR 
THE SYSTEM WITH LIMITING (FIGURE 7) 


Path Tangent Curves 


As previously mentioned, only the discrete 
points in the A-plane which represent the 
sampled quantities are significant, and not the 
points in between. Ifthe consecutive sample points 
are connected with line segments, the resulting 
polygonal line will be referred to as a solution 
path. The slope of the path segment is Ay/Ax 
where y = Ax. 


A method for tracing solution paths for equa- 
tions which may have continuous nonlinearities is 
obtained by considering the similar differential 
equation. The ordinary phase plane trajectories 
ofthis similar equation are tangent to the segments 
of the solution path. The similar equation is ob- 
tained by replacing Ax by x and A¢x by x. The 
slope of the ordinary phase plane trajectory is 


dy/dx and the slope of the path segment is Ay/Ax. 


The usual methods of obtaining the phase plane 
trajectories may thus be used to study the sampled 
system. For the previous example, the similar 
equations are: 


ous 4d Jc + c| >1 
ee ; (29) 
Co bace +. CG = 0 le + c| <a 

The path tangent curves are sketched in 
Figure 9 with a typical path indicated. 


A path is obtained by continuing along the 
tangent of the path tangent curve through the 
starting point until the x-coordinate has changed 
by the Ax indicated at the starting point. This 
procedure becomes indeterminate for starting 
points where x = 0. Along the x-axis, the original 
equation is used to determine Ay. 


For systems which switch on some boundary 
in the A-plane, the path tangent curves switch on 
the same boundary and do not overlap as the tra- 
jectories may do. The singular points of the 
similar differential equation are also singularities 
in the A-plane. Perturbation about these singu- 
larities gives information about their stability. 


The system shown in Figure 10 is described 
by the difference equation 


4 


Sa 


FIGURE 9.PATH TANGENT SOLUTION 


Peesiee a ihhitees 5 SO (30) 
or 
2 2 2, 
Nex. tee xmas Ax (2 scmaidiisc 9 = 0 (31) 
The similar differential equation is 
ase (2 . 2 
x + (2+x -1)x + (2+x -1)x = 0 (62) 


FIGURE 10.SYSTEM DISCRIBED BY 
Aex + (24x27 -1)ax + (24+x°-1) x=0 


Using the isocline method, the path tangent 
curves are obtained as shown in Figure 11. 


Several paths are indicated. Path A is divergent; _ 


path B results in a periodic solution. 


Conclusions 


We have emphasized here the use of differ- 
ence equations to describe sampled systems. 
From this emphasis, the A-plane as a means of 
analyzing sampled system behavior evolved. 
Much about the method remains to be developed 
(for example, the possibility of its use to predict 
limit cycle behavior). As it stands, though, the 
method furnishes further insight into the behavior 
of sampled systems, and in particular, those 
containing nonlinearities. 
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ON THE EXISTENCE AND UNIQUENESS OF THE 
OPTIMAL MULTIVARIABLE SYSTEM SYNTHESIS 


Mihajlo Mesarovié 
Case Institute of Technology 
Cleveland, Ohio 


Many problems in the study of multivariable 
systems are not of the same nature as the prob- 
lems normally associated with single variable 
systems, This has interesting consequences 
when solving practical problems of the synthesis. 
For instance, when in the synthesis of single 
variable systems the available analytical pro- 
cedures cannot provide the solution one can 
simulate the problem on the computer (analog or 
digital) and obtain numerical solution for the 
particular situation at hand. However, in many 
problems in synthesis of the mltivariable 
systems the solution itself depends upon the 
way of representing the system. Such a situatim 
appears in synthesis of the cross controllers 
when the resulting changes in the interrelations 
of the closed loop system should be taken into 
account. 


We shall consider some specific properties 
of the multivariable system synthesis which are 
valid for a broad class of problems rather than 
look for the closed-form solution of some 
special problems, Actually in the absence of 
an adequate closed form solution one can look 
for a numerical solution but in doing so one 
should be aware of the general properties of 
the problem at hand. In a monograph to appear 
soon one can find more detailed discussion on 
the specific properties of the multivariable 
systems.~ Here we shall be concerned with the 
problem of the feasibility of synthesis. More 
accurately we shall consider the following 
problem: Given a set of input variables and a 
set of the corresponding desired output variables 
which kind of systems might offer a solution and 
whether the resulting system is the only one to 
satisfy the requirement of the synthesis. 


In the conventional formulation of the 
synthesis problems in single variable systems, 
the existence and uniqueness of the solution 
ultimately depends upon the possibility of find- 
ing a unique solution of the equations defining 
the optimal conditions. However, if the 
synthesizing procedure does not offer a unique 
solution it is important to find the reason for 
that. In general, there are three steps in 
which the optimizing procedure can fail: 


1) The equation for optimal behavior may 
have no solution, i.e. even the ideal optimal 
system cannot be found. 


2) The ideal optimum can be found but not 
a satisfactory approximation of the optimum in 
a restrictive class of functions such as a class 
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of physically realizable functions technically 
realizable functions (in a given sense), etc. 


3) One cannot synthesize the system with 
the hardware at hand, even if its behavior 
belongs to a restrictive class of functions. 


The problems in each of the above three 
steps become even more complex in the case of a 
multivariable system. In addition there is one 
step more in which the unique synthesis of the 
multivariable systems can fail: namely, the 
optimal conditions can fail to offer a large 
enough set of equations, the solution of 
which should specify the optimal system. This 
step precedes the above mentioned, and in what 
follows we shall be discussing relating problems. 


Binary Representation of Systems 


Before further discussion let us consider 
the possibilities for a more precise definition 
of the synthesis task. In the synthesis of the 
system mainly two kinds of problem may appear: 
(1) The inputs and desired outputs are given and 
one should synthesize the system, This class of 
problem could be termed filter-problem. (2) In 
addition to the desired outputs and inputs the 
system (or process) is given and one should 
synthesize another set of inputs, often termed 
manipulated variables in the control type of 
problems. 


We shall show in this section in what 
sense these two kinds of problem are analogous. 
To this end we shall use binary representation 
of the system. In this kind of representation 
every function among the variables is repre- 
sented as a binary operation, i.e. it depends 
upon two variables (two inputs) and produces one 
variable (one output). The linear function can 
now be represented by convolution as binary 
operation which, as two inputs, has the actual 
input and weighting function of the element. A 
static nonlinear element can be represented in 
many different ways, It can be represented by 
multiplication as binary operation with one of 
the multiplicants (which is now time variant) 
determining the kind of nonlinearity. Multi- 
plication is not, in fact, always convenient 
because it occasionally requires infinite 
values for one of the multiplicants, As a 
better way of representing the nonlinearity the 
potential with time variant exponent can be used, 
The memory-type nonlinearities might require many 
cascading binary operations for an adequate 
representation, 


The choice of the kind of binary operation 
is to some extent arbitrary. Preferably 
the operation itself should not include quanti- 
tative information, or at least no information 
subjected to synthesis. In other words the 
binary operation should specify only the kind of 
operation to be performed while all associated 
data and variables should be represented as 
inputs. A linear system e.g. can be represented 
by the binary structure of the system which 
includes only convolution and addition as binary 
operations and two sets of inputs. One set of 
the inputs consists of weighting functions of 
all elements. They could be regarded as gener- 
alized parameters and will be referred to as 
behavioral inputs, The other set consists of 
the input variables of the system and for the 
sake of distinction they will be termed environ- 
mental inputs. A general representation of the 
system is as in Fig. (1). The binary structure 
includes many different binary operations. 


The problem in the synthesis that we want 
to consider can now be formulated in the follow- 
ing way: Given the binary structure of the 
system and one set of inputs (environmental as 
well as behavioral inputs can be included in 
this set) the remaining set of inputs should be 
synthesized. 


A General Form of Necessary Conditions For 
Optimal Multivariable Synthesis 


To discuss problems in the first step of 
synthesis we are in need for necessary conditims 
under rather general circumstances, To this end 
we shall take as the performance criterion, the 
functional of the form 


t 


CecoghA * FACE oSese 7,(*) | at 


by 


(1) 


where: Fis a nonlinear function of all the 
outputs Yo sores Vys 


The system is assumed to be described by 
the equations 


Gs [ays By, Bley), one tes By Bly)s ov 


(2) 
oce 4? jy, P(y,)...-- Ins Sins By) | = 0 


where Z ++ 2% are inputs to the binary 
structure which should be synthesized B(y4) nee 
tee ay) are functionals which describe the 


linear parts of the system. The inputs to the 
system X19 eee X, are supposed to be known 


deterministic time functions and therefore are 
not explicitly shown in &q. (2). 


Consider the case where the optimization 
should be performed only with respect to one of 
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the inputs, say,z,.. The calculus of variation 
technique gives fér the optimizing condition 


ZA * G.(z,) = 0 (3) 


k 


where the operator T defined by 
OK(x) _d_ /dK(x) OK(x) (of) _ 
oie a erm eg rer, per gh 


ox 
was introduced. 
either G or F, 


Ak 
m 


(4) 


Function K in Eq. (hk) can be 


where : A ic 


Determinants A and A do not depend upon the 
input z,. which is to be synthesized, 


Jj 
oS. {Ty i @.(y,.)} 
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Another interesting form for the necessary 
condition of the optimum can be obtained if the 
behavior of the system is given by the equa- 
tions. 


(5) 


be hy lz, a? as) = “ARE a? 


Ma, Ds Jy i> By)» eee ¥5-1? I5-1 


Bly;_4)s Vj+? Tye B(Y,,) CAO n? a 
Ay.) ] (6) 
The optimal condition then becomes 
oy. oy 
oF 4] j 
i 35 a = 0 (7) 


Proceed now to consider the case with the 
random inputs, Assume that the system is linear 
and adopt the cononical structure for it. The 
equatiors describing the behavior of the process 


are 
cs 0) 


¥,(t) = z i Psy () X(t - +) dv (8) 


Take as a performance measure a function of 
all the out puts Y,, eae ty and reference 


variables Y. 


lr? eee Yor 


Y (9) 


e=e i coe ty ww °* Tee 


Expanding the performance criterion in 
Taylor series and taking only the first members 
one obtains 


e= e(y,5 e'eie Yn? Vyps eee Yeo) + 
ntnr cS) e(y,, eee Vn? Yr? coe Yor) 


tes 
1 ) yy 


ES + x] 


where bar denotes mathematical expectation of 
the respective variable. 


(10) 


The mean square of the performance variable 
is now 


ntnr ntnr 
2 de de 
eo= 5 2 xe xe G., (0) (11) 
j x OY; % IK 


where Bs. is the cross-correlation function. 


As given by Eq. (11) the final performance 
measure is a linear function of various correla- 
tion function of the output variables and 
therefore has no extremum in respect to them, 
For greater gengrality assume that the perfor- 
mance measure e“ is a nonlinear function of 
correlation functions 


s?er[... By (0) oe | 


The system can now be represented by a 
deterministic model which has equivalent behavior 
when measured by the performance criterion (12). 
The equivalent deterministic model has hn 
outputs and n@m? inputs, The problem of 
synthesis is now in determining the behavioral 
inputs, i.e. the transfer functions of the 
elements. 


(12) 


Existence and Uniqueness Numbers 
of Multivariable Synthesis 


We shall now discuss the feasibility in 
formation of optimal conditions, As was 
mentimedin section 1 this is the first step in 
the multivariable synthesis in which either 
existence or uniqueness of the synthesis can 
fail to be obtained. 


The ultimate goal of the synthesizing 
problem, as defined here, is to determine a set 
of inputs to the system. However to achieve a 
total optimum according to Eq. (1) only one 
subset of the inputs could be sufficient. The 
inputs belonging to this subset will be called 
the bounded inputs because they are completely 
specified with the synthesizing procedure, The 
remaining inputs of the set will be termed free 
because no matter how one chooses them the 
bounded inputs might assure a totally optimal 
system. 


The existenee and uniqueness in formulation 
of the optimal conditions can now be determined 
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in an extremely simple way by considering only 
two numbers. The first one, termed existence 
number, represents the difference between the 
number of the outputs and the number of the 
bounded inputs. 

(13) 


D=r q7 0 
Its sign indicates the existence of the 
set of the conditions which specify the optimal 
system. The second number, termed the uniqueness 
number, represents the difference between the 
total number of the elements in the input set and 
the number of the bounded inputs or, in other 
words, it represents the number of elements in 
the set of the free inputs. 


(1s) 


The usefulness of the concepts of existence 
and of uniqueness number comes from the fact that 
in a given situation they reflect the kind of the 
systems structure and at the same time the 
necessary information one should have for the 
optimal synthesis. 


Eats =" 54 


For the linear systems one has the condi- 


tions (7). Taking into account that for linear 
systems one has 
oy. 
Jp 


the existence and uniqueness numbers become . 


D= 0;R=n(m-1)>0 (16) 

Therefore the conditions for the total 
optimum of the linear system exist but are not 
unique, i.e. many different conditions may be 
formulated resulting in many different systems 
all having some behavior (as measured by the 
given performance criterion). 


For nonlinear systems,one obtains from 

Eqe (3). 
D=0 R>O (17) 

The optimal conditions exist and might 
also be unique. The nonlinearity in the 
structure of the system can be measured just by 
R, Whenever one has R> 0, i.e, there is at 
least one of the inputs free, the system is not 
completely nonlinear in structure, The degree 
of the nonlinearity in structure could be 
measured by the value oP R, 


When the given inputs to the system are of 
random nature, one obtains according to Fig. (2) 


(18) 


D is now negative and the set of conditions 
for a totally optimal system does not exist, 


D = n(m - 3n) 


This comes from the fact that a deterministic 
equivalent model has more outputs than inputs 
available for the synthesis, The existence of 
the optimal synthesis depends now upon the 
correlation of various outputs and reference 
variables, e.g. if every output is correlated 
only with the corresponding reference variable 
one has 


D = n(m = 2) (19) 


The existence of the optimal synthesis 
depends now upon the number of the variables 
associated with the system. One can obtain the 
optimal system only if n= 2, The increase in 
the number of variables of the system enables 
an easier synthesis, 


R and D reflect the basic properties of 
the system in an utmost simple way and therefore 
they should be considered as characteristics of 
the system as e.g. the numbers of inputs and 
outputs. Therefore for the synthesis of a 
maltivariable system one should have in 
advance: 


1) Four characteristic numbers of the 
system, n, m, R and D, 


2) A set of inputs and a set of the 
respective desired outputs. 


3) The specification of the binary 
structure is perferable and in many practical 
cases presents no limitations whatsoever, 


Additional Specifications for Synthesis 
of Totally Optimal Systems 


For systems, R>0, additional specifications 
are necessary if the synthesis is to be unique. 
This is in sharp contrast with the properties 
of the single variable problems. 


Two ways might be offered for specification 
of the additional conditions, 


1) One could specify in more detail the 
structure of the system as e.g. by requiring 
a given kind of interactions inside the system, 


2) One can specify a large enough number 
of sets of inputs and desired outputs. One has 
the possibility then to cover an entire class 
of input functions offering a more realistic 
synthesis of the system, 


The synthesis with the additional specifi- 
cations presents a basic and very interesting 
problem in the multivariable study but is out 
of the scope of this presentation. 


lu, Mesarovié » "Foundations for a Control 
Theory of Multivariable Systems", Technology 
Press, John Wiley, 1960. 


Fig. 1. Binary representation of a multivariable system. BS fy 
binary structure, Z; = environmental inputs, Xj = behavioral 
inputs, e = performance measure. 
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ON OPTIMAL AND SUBOPTIMAL POLICIES IN THE 


CHOICE OF CONTROL FORCES FOR FINAL-VALUE SYSTEMS 


by 


Masanao Aoki 
Department of Engineering 
University of California, Los Angeles 


Summary 


A simple control system whose state var- 
iable xy is described by the difference equation 
(1) is considered. 


Kya 2 Sy iateh vn. os tear 


6<a <li 


where vy is the control force of the system and 
pn is the Bernoulli random noise with the prob- 
ability parameter p. 


If $ (xy) is the performance index of the 
final-value system, then the problem is to find 
a sequence of control forces {vn Peak Pca bead A ea Nb 
which minimizes the expected value of ¢ (xy). 
An optimal sequence of v, is determined by sol- 
ving a recurrence relation on n of the criterion 
function h,, (x;p) defined as follows: 


Min E [6 (xp)] 


7+ Vy 


hy (x;p) = Min Min... 


Vi V2 
= the minimum of the expected value of 

¢ (xn) obtainable in the n-stage control 
process by employing an optimal policy, 
starting from the initial state variable x 
when the probability parameter of the ran- 
dom disturbance is p. 


The recurrence relation is obtained by the 
usual application of the principle of optimality 
of dynamic programming technique. 


It is proved that hy (x;p) = hy (-x; 1-p) 
holds and this fact is used in reducing by half 
the amount of computation when it is necessary 
to solve the recurrence relation numerically. 


If the values of control variable v,, (x;p) 
are restricted to 1 and -1 as in contactor servo 
systems, the boundaries between 1 and -1 con- 
trol forces become too complicated to be deter- 
mined analytically except for few special cases. 


By solving vy (x;p) computationally for the 
¢(xyx) = xn? case, it is seen that vy (x;p) agrees 
with v, (x;p) fairly well for alln> 1. 
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v1 (x;p) is the optimal control force when there 
remains only one chance of exerting control for- 
ces. Hence, the suboptimal policy of always ap- 
plying control forces as if only one more error 
correction is possible may be expected to be 
close to the optimal policy in minimizing E(x’). 
The control force v, (x;p) is linear in x and p and 
is a much more simple function to mechanize. 
This suboptimal policy is tried by the Monte Car- 
lo method and found to be only slightly inferior 
to the optimal policy. 


The behavior of xn is investigated by assum- 
ing the control forces are given by the subopti- 
mal policy. This approximate analysis should 
be good in view of the agreements between the 
optimal and suboptimal policies. 


It is important to realize that if the adop- 
tion of a suboptimal policy results in a simpli- 
fied mechanization of the optimal control forces, 
with only a slight reduction in the system perfor- 
mance, then the suboptimal policy might be opti- 
mal in a certain enlarged performance index. 


What seems to be the most desirable ap- 
proach to engineering problems is a unified or 
well-integrated one where the analytical and 
computational algorithms are used to supplement 
each other. In this paper, an attempt is made to 
illustrate the point, presenting at the same time 
a new approach to analysis and synthesis of a 
certain class of control systems. 


introduction 


Time behavior of a control system can be 
regarded as successive transformations on a 
vector, called the state vector, which represents 
the state of the control system as a vector in 
phase space. 122238 


Types of transformations to be applied are 
determined by specifying the control vectors as 
functions of the present state vector and time. 
Thus, if the time is assumed to advance in dis- 
crete steps, a problem of determining proper 
sequences of control vectors, given the criterion 


of performance, can be formulated as a problem 
of multi-stage decision processes. The deci- 
sion at each stage is a proper choice of the con- 
trol vector from the domain of control vectors 
which may or may not be known in advance. 

That domain is, however, assumed not to change 
with time. 


Functional equation techniques of dynamic 
programming are well suited to treatment of 
control processes as multi-stage decision proc- 
esses. 


By solving a functional equation obtained 
from the application of the principle of optimal- 
ity to the control process in question, a sequence 
of optimal control vectors is obtained. This 
sequence is called an optimal policy. Except in 
a few isolated cases, analytical solutions of 
functional equations are generally not available. 
Thus, in order to obtain an optimal policy, it is 
often quite necessary to resort to some compu- 
tational techniques to solve the functional equa- 
tion numerically or to solve the functional equa- 
tion in some approximate iterative fashion and 
make the sequence of approximate solutions con- 
verge to the true solution. 


Even in numerical solutions of a functional 
equation, because of the limited capabilities of 
existing digital computers, it may be feasible to 
solve the equation only approximately. Thus de- 
pending on the complexities of functional equa- 
tions, it often becomes necessary to investigate 
a sequence of non-optimal decisions. It is called 
a suboptimal policy. 


Even when an optimal policy is available, 
we may be interested in some policy slightly 
inferior to the optimal policy, if the adoption of 
this suboptimal policy results in vastly simpler 
mechanization of control devices in a control 
system. 


In the following, this concept of optimal 
policy versus suboptimal policy is pursued by 
describing the situtation where a stochastic con- 
trol process serves as an approximate solution 
to a certain adaptive process and further by in- 
vestigating a certain simpler control policy in 
the stochastic process which serves as a sub- 
optimal policy. Concepts of stochastic and adap- 
tive processes will be made precise later. 
Roughly, in stochastic processes, we are given 
distribution functions of random variables, 
whereas in adaptive processes, relevant distribu- 
tion functions are imperfectly known, making it 
necessary to gather additional information on 
distribution functions from the history of the ob- 
servations on random variables. 


Adaptive Control Processes 


Distribution functions for random vari- 
ables affecting the performance of time~dis- 
crete control processes are assumed to be mem- 
bers of a certain parameter space, which is as- 


W72 


sumed to be discrete and finite. A random vari- 
able R, therefore, has the distribution function 
of the form F (r, @) where 


p (@= aj) = 2. 220, red Redeye A ie 
k 
» Zoe 
ileal 


Let us associate the states of nature with each 
aj, i=1, 2...K. They are assumed not to 
change with time. 


If zi = 1 for some i, then the control system 
is known to be in the ith state of nature. This is 
a stochastic situation since the distribution func- 
tion of the random variable R is given by F(r, aj), 
a known function, 


If,0 <2; < 1 for all j--.1.¢2 .5.K,) ten ghe 
control system is adaptive. Depending on the 
history of observations on R, the estimates zj 


are assumed to be modified by the Bayes rule. 


The adaptive model described above is a 
special case of the adaptive models of Bellman 
and Kalaba. 4 


The difference equation of a final value con- 
trol system S is given by 


xnt1 = T( Xn, rn, vn), n=0, 1...N-1 (1) 
where xy is the state vector and vy is the control 
vector of S. Let &xy) be the criterion of this 
final value control system. 


The random variable ry is assumed to be 
independently and identically distributed with the 
distribution function F (r, a). 


The control system S must decide to apply 
vn} which minimizes the estimated expected 
value of $(xy) , given the present state vector 
x, the number of the remaining decision stages 
n and the present estimate of the probability 
P( aie of) ata pel a 25.2 Kronen — (2472 Zi). 
The estimated expected value of $(xy) employ- 
ing the optimal sequence of v, is expressed by 
ky (x;z) where x, z and n are as described above. 


In the following, for the sake of brevity, we 
will restrict ourselves to the case of K = 2. 


Let us redefine z as 


Pe sayy = 2) P le sa, ) 21-2 

Then, if we define hn (x;aj), i= 1, 2 cor- 
respondingly to kn (x;z), as the expected value 
of ¢ (xy) following an optimal policy, then ky (x;1) 
and k,, (x;0) will become equal to hp(x;a;) and 
and hy (x;@,), since once z = 1 or z = 0 is real- 
ized, the value of z will not change. 


To obtain an optimal policy, it is necessary 
to solve a functional equation for kn (x;z). Under 


the above conditions, however, it can be proved 
that ky (x;z) is concave in z, 


Ky (x3z) = z kn (x31) + (1-z) ky (x30) 


= Z hn (x3@, )+ (1-z) by (x;e,) (2) 

This equation means that knowledge of 
hn(x;aj), i= 1, 2 provides a lower bound on 
ky(x;z) without solving the functional equation 
for kn(x;z). 


Noting the fact that aj in hp (x;aj) is a pa- 
rameter and not a part of the independent vari- 
ables of the functional equation as z is in ky(x;z), 
it is seen that the functional equation for hy(x;aj) 
will in general be easier to solve. 3 


The knowledge of an optimal policy and the 
value of ky (x;z) is equivalent, i.e., knowing 

Kn (x;z)} , all optimal policies are obtainable. 

n the other hand, the knowledge of an optimal 
policy suffices to generate {kn (x;z)}. Therefore, 
the approximate knowledge on ky (x;z) as sup- 
plied by (2) may be used to provide a suboptimal 
policy which is constructed from our knowledge 
of the corresponding stochastic control process 
alone. 


The usefulness of (2) depends entirely on 
how good this approximation is. We will come 
back to (2) later when the numerical results are 
discussed. 


Stochastic Control Processes 


In order to illustrate several points made 
in the previous section, let us specify the under- 
lying difference equation (1) for S further to that 
given by (3). 


Roe ene tat Va Oe 2 (3) 
where ry is taken to be the Bernoulli random 


variable 


with probability p 
with probability 1-p (4) 


The criterion of performance is taken to 
be 
-1 
E (Jy) = E (xy? + A . Wa") (5) 
i=o 


where the 2nd term is included to express the 
situation where there is a constraint iene BAe total 


amount of control force in the form ye vy" T< V. 


10) 


Then, defining h,, (x;p) to be the minimum 
of the expected value of Jy obtainable in the n- 
stage control process by employing an optimal 
policy, starting from the initial state variable x 
when the probability parameter of the random 
disturbance is p, the recurrence relation for 
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hn (x;p) is given by 


h, (x;p) = Min [A Vo? + p ( xae)? + (1-p) (x_)| 
Vo 


hn+i (x;p)= Min [a Vn? + P hy (x45p) 


3 +(1-p) hy (x_sp)] 


Deh lee ere Nicol! 


(6) 


where Roe 


3 eo) 36 ae As) se Ep 


a_i 20 96g) Joy ce ie) 

If no other constraint is imposed on v, then, 
as expected, h, (x;p) is quadratic in x and v,(x;p) 
isslinear in, x, 

hy (xsp)= (b? - B?) 2 tA Nala eae 

At+1 
eae (n-1) 

NEL) of +. J yelt-2) 


‘ > 
A+ 1+ a? +2 (n-1) 


ciiele se 


PE te 


[oP x+B(ltat..te 


an-lienxy +B (1 +a¢...¢ an7l) 


Feb) se yo (8-1) 


Ne Ieee hee. 


where b = E (b) = (2 p- 1)b. 


However, if vy is constrained to be 
(8) 


as in a contactor servo system, then explicit ex- 
pressions for hy (x;p) and vp (x;p) are no longer 


Vn = Mor-m, fag) 224s) 


available. Since Y nee 
1=0 
performance can be taken simply to be 


E (x17?) 


Nm2, the criterion of 


The recurrence relation is now given by 
hy (x;p)= Min [p (22a) (x_)?] 
Vo= +m 


hn+1 (x;p)= Min [ phn (x,;p) + 
Vn= +m 


(1-p) by (x_;p)] (9) 


n= ils 2 eee N-1 

Although explicit expressions for hy and vn 
are not available, it can be shown by inductive 
arguments that hy (x;p) = hp (- x;1-p) holds.’ 


Due to this symmetry, the amount of com- 
putation necessary to solve hp (x;p) for a given 
range of x and 0 < p< 1 is reduced by half. 


The symmetry property is seen to hold 
even when the domain of vy is expanded to in- 
clude more than two values so long as for vp = a, 
the corresponding value vy = -a is in the domain. 
(vn = 0 is the special case of this). 


If the value of p is regarded as the para- 
meter a of the previous section, i.e., 
P (p=pi]=z, —-P. [p=p,]= 1-z 
il > Pi > P,2 0 


then ky (x;z) as defined in that section will satis- 
fy the recurrence relation 


k, (x, z) =Min 
Vict ae gay 


[Pr o(x4) +(1-ps 6(x_) 
+ (1-z) [P. @ (x4) + 
(1-p,) ¢ (x_)| 
Zz [ps kn (x4, z') + 
(1-p,) kn (x_, z" |+ 
(1-z)[p, BAAxe, 2')i+ 
(=p) cn: (> | 


kn+1 (x1 z)=Min 
Ve=tm 


where x,and x_ are as defined in (6), and 
Pi Z (1-p;) z 
Pi z+p,(1-z) 4a (1-p; )z+(1-p,)(1-z) 


The recurrence relations (9) and (10) are 
solved computationally with the boundary condi- 
tion that |x,| < Dfor all n for the same system 
parameters. This boundary condition means 
that there is a physical stop beyond which the 
system cannot go. 


ites 


FE 


By comparing ky (x;z) and z hy (x;p,) +(1-z) 
hn (x;p,) the approximation of (2) is tested. 
Figure (1) shows kp (x, z) as a function of z for 
several n and x, 


Since the approximation of (2) amounts to 
regarding ky (x, z) to be linear in z, the degree 
of approximation is seen to be better for larger 
aay 


One-Stage Suboptimal Policy 


It is shown in previous sections that an 
optimal policy for the stochastic system (3) is 
obtainable by solving (9), with knowledge of 
{hn (x;pj)} i = 1, 2 presupposed in the approxi- 
mation of ky (x, z) in terms of hp (x;p;,) and 
ice (x;p,). 


Let us next turn to the investigation of the 
approximate solution of (9), i.e., the determin- 
ation of a suboptimal policy of the stochastic 
control system (3) itself. 


Figure (2) shows an optimal control vari- 
able as a function of x and n for p=0.625, D=4, 
a= 7/8, b= 1/16, m = 9/128, n < N = 12.* Not- 
ice that the boundaries between vp (x;p) = +m 
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and vn (x;p) = -m is not a simple straight line. 
When N = 1, the optimal control vector is 
given by 
v; (x;p) = - m:sgn ( x + (2p-1) b) (is) 


The switching boundary, then, is a straight 
line given by 
- (2p-1)b/a@ (12) 


x= 
and shown in Figure (2). 


If (11) is used as vn (x;p) for n> 2, then it 
constitutes a suboptimal policy. 


From Figure (2) it is noted that for x > 0 the 
optimal and suboptimal policies agree fairly well 
although the agreement is not particularly good 
forests 


It might be expected, therefore, that if the 
system starts from the initial position x = D, 
then the suboptimal is a fairly good approxima- 
tion for the optimal policy. 


To see to what extent the conjecture is cor- 
rect, the Monte Carlo method was used to follow 
the system behavior from initial positions xo 
+ D for 40 trials each. 


With x, = D, E (x?) optimal = 0. 00367, 
E (xq) suboptimal = 0. 00374. 
With x, = -D 
E (xy?) optimal = 0. 00437, 
E (x}q2) suboptimal = 0. 00441. 


Despite the relatively similar values in 
E (xy?) , the suboptimal policy for x9=D appears 
to be better as conjectured, by the fact that out 
of 40 trials with x, = D, the optimal and sub- 
optimal policies gave the same xy? in 21 trials, 
but in similar 40 trials with x, = -D, the opti- 
mal and suboptimal policies did not give the 
same xy? in any case. 


Although these results are by no means con- 
clusive, they tend to support the conjecture that 
the one-stage suboptimal policy is a fairly good 
approximation to the optimal policy for the con- 
trol system under consideration. 


The adoption of this suboptimal policy simp- 
lifies the control problem considerably. 


The approximate analysis of the control 
system with the one-stage policy will next be 


* Nis taken tobe 12. For N> 12, the system 
stays pinned to either x = D or x = -D because of 
the particular boundary condition of the system 
until the number of remaining stages, n, be- 
comes less than 12 or 11. 


presented. 


The Control System Behavior 
With One-Stage Policy 


The state variable under the suboptimal 
policy is given by 


Xnti = @Xn+ n+ Vn (3') 
where vn = -m‘sgn (a Xn + (2p-1)b ) 
Ry { _+b_ with prob. p 
=b» with prob. 1=p 
USS lols seal SE 


Let us note that if x, > 0 then xnii < Xy 
and if xn < xc, then xnii1 > Xn where xc (2p-1). 
b/a, because (m-b)> 0 by assumption. Namely, 
if the system deviation from the origin ( the 
position of the equilibrium in the absence of the 
random disturbance r) is positive, the deviation 
at the next stage will invariably be smaller than 
the present one. If the present deviation is less 
than xX,, then the system will move towards the 
origin at the next stage. 


For Xe < ¥n <0, Vn =-m 


ADU tase Keb 10 - (m + 2pb). 

That is to say, once Xj becomes greater 
than xy, = Min[- (m+ 2pb), - (2p-1) b/a@], then 
Xx, cannot become smaller than x,,, for all k>n, 
i.e., xp, is the lower bound. Similarly, one 
sees that 


Xu =m + 2b(1-p) 


provides an upper bound. That is to say, if xn 
< x,, then x, < xy for all k>n. 


Figure (3) shows the typical system behav- 
ior using the one-stage suboptimal policy. 


One noticeable feature of the system be- 
havior is that xp approaches the origin for a few 
stages consecutively and then xn jumps away 
from the origin and from there it again ap- 
proaches the origin for the next few stages. 


Let us now investigate the number of 
stages before a jump occurs. This serves to 
show how the jerkiness of the system depends 
on system parameters, a, b, m and p. 


Suppose that xj, = -(m + 2pb), and Xo = Xj,. 


Let k be the first time that x, > X¢ is 
realized. 


The largest k will be realized when ry is 
equal to - b for all k stages and the smallest k 
will be given when rp = + b occurs consecutively. 


In the former 
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(m-b)(1-ak) (2p-1)b 
Xk F ~ ak x0 + SPS eee ee (Gis) 
In the latter 
eon (m+b)(1-aK) _ _ (2p-1)b 
X; ak xo + aa aa = es (14) 
Pecoman( 143) 
m-b 


in [= + eee) - In[m+2pb+ Taal 
In [0] Hoel Magee (15) 


k< 


From (14) 


b 
dnl. 
ree fee 


(2p-1)b] © m+b 
aP i | In[m+2pb He Riles 
n (2) 
With a = 7/8, p = 0.65, b = 1/16 and m = 9/128, 
k=5 in (15) and k=1 in (16). k = 5is realized 


with probability (1-p)5 and k=1is realized with 
probability p. 


| 
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If x9 = Xy, then the similar argument gives 


rey ee : @p- Up) In[m+2p(1-p)+722I 


A {] 
(17) 


the maximum k is given with - sign and the mini- 
mum k is given with + sign in (17). The maxi- 
mum k occurs with prob. pk and the minimum k 
with prob. (1-p). 


Thus, after n becomes such that xy, < Xp < 
Xy is realized, then the bounds on the number of 
consecutive stages, k, where the system ap- 
proaches the origin monotomically are obtained. 


Conclusions 


It was shown that for a class of final value 
systems considered in this paper, ky (x, z) is 
concave in the parameter z which represents the 
a priori knowledge of the states of nature. 


Policies derived from the functional values 
z hn (x;p,) + (1-z) hy (x;p,) are shown to serve as 
initial suboptimal policy for adaptive control 
systems. 


When the control variable is of the binary 
form +m, no explicit dependence of hp (x;pj) and 
Vn (x, pi) on x is obtained. 


However, one-stage suboptimal policy was 
shown to be a fairly good approximation to the 
optimal policy and was used to establish lower 
and upper bounds on the system deviations. 
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priori probability, ¢. 


176 


THE NUMBER OF REMAINING 


CONTROL STAGES 


OW MONO VU PWN 


pec 


THE STATE VARIABLE OF THE SYSTEM, X -® 


d 
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Fig. 2. Boundaries between +m control force and -m control force in the optimal and the sub-optimal 


policies at p=0.625. 
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A STUDY OF ASYNCHRONOUSLY EXCITED OSCILLATIONS 
IN NONLINEAR CONTROL SYSTEMS 
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Summary 


This study was initiated as an effort to 
explain certain oscillatory phenomena observed 
in an aircraft and weapon control system 
operated in a particular mode. The oscillations 
in question were of the undamped limit-cycle 
type and their presence could very clearly be 
correlated with the degree of noise corruption 
of the signals. 


It is demonstrated in this paper that cer- 
tain types of nonlinear systems although being 
inherently stable may be driven into oscilla- 
tory modes not only by random signals but also 
by any high frequency periodic or nonperiodic 
signal possessing a certain energy content. 
is of interest to note that the tendency for 
hunting and also the hunt frequency usually are 
completely independent of the frequency of the 
excitation signal, i.e., the phenomenon is of 
asynchronous nature. 


It 


A general theory, the validity of which 
has been tested by both analog and digital means, 
is presented and utilized to demonstrate how 
this phenomenon may be predicted from informa- 
tion on circuit data. 


Introduction 


Although the phenomena to be discussed in 
this paper can only occur in systems incorpo- 
rating certain types of nonlinear components, 
it should nevertheless prove enlightening to 
take a quick look at some of the characteristic 
features of linear systems. This will be done 
in order to demonstrate the uniqueness of the 
phenomenon and also in the hope that it may be 
possible to employ some of the well-established 
linear analysis techniques in explaining it. 


A system is defined as linear if its 
dynamical behavior can be described by a 
system of linear differential equations, par- 
tial or ordinary. From the analyst's point of 
view, a linear system is quite attractive be- 
cause of the relative ease with which its 
behavior can be predicted. The control engineer 
in particular is appreciative of the following 
two characteristics: 


ek The stability of a linear system is 
uniquely determined from the parameters of the 
system only. 
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2. In response studies the super-position 
rule can be employed. 


These two features, probably to a higher 
degree than any others, put the linear system In 
sharp contrast with the nonlinear one. The first 
one of the above rules states that a linear 
system cannot be driven into periodic or growing 
oscillations by any combination of external in- 
put signals of finite duration. One may, in 
contrast, quite simply demonstrate a multitude 
of nonlinear systems which, for Instance, are 
stable for low signal levels but go into limit 
cycles (hunting) should they be excited by high 
level inputs. 


The second rule points out a very important 
"'nonmixing!! characteristic of a linear (stable) 
system which indeed constitutes the basis for 
all the well-known frequency response techniques 
in use. A linear system will at its output ter- 
minals contain only those frequencies it has 
received at its input. A nonlinear system, on 
the other hand, never displays such a frequency 
discrimination. It may accept a single or a 
few frequencies at its input and after a compli- 
cated mixing process deliver an infinite, but 
discrete, spectrum at [ts output. The new fre- 
quencies created In this process are referred to 
as a modulation product. J|t will prove advanta- 
geous in the following to consider the nonlinear 
components as modulators. 


It has been known for considerable time 
that injection into a control system of a high 
frequency dither signal will in certain instances 
upgrade the performance of the system. The fre- 
quency of this signal is usually chosen consi- 
derably above the cutoff of the system so that 
its presence cannot be detected in the output. 
The "dynamical lubrication'' provided by such a 
signal was first pointed out by MacColl’ who 
analyzed its effect in an on-off type servo. It 
Is today quite common to diminish the effect of 
static friction by use of a superimposed high 
frequency signal. (This ‘lubrication effect" 
is explained in a simple fashion by the theory 
later to be developed in this report.) 


MacColl also pointed out how, by Inserting 
a high frequency signal into the on-off servo 
mentioned above, one could improve its sensiti- 
vity for small signals. Oldenburger had dis- 
covered experimentally that a high frequency 
signal superimposed upon the Input signal of a 


nonlinear governor control system tended to 
improve its stability. Indeed it was possible 
to completely quench self-oscillations under 
certain conditions. He reported his findings 
in two separate papers2»3, in the latest of 
which he attempts to analyze a particular single- 
loop system with a well-defined nonlinearity in 
the forward loop. This method of stabilizing 
nonlinear systems subject to limit cycle oscil- 
lations is known under the name of ''Signal 
Stabilization". 


The topic of the present investigation is 
to study the opposite problem: Is it possible 
to excite an inherently stable system by a high 
frequency noise Signal in such a manner that 
the system ceases to be stable? To the author's 
knowledge, this problem has never concerned con- 
trol people. As a matter of fact, an extensive 
search of the literature on nonlinear phenomena 
reveals very little work in this area. It 
should be made clear at this point that the 
phenomena discussed here (‘''Signal Stabilization! 
included) differ essentially from the well- 
known phenomenon of so-called nonlinear or 
"subharmonic resonance' in which the existence 
of a rational ratio between system frequency and 
frequency of excitation always must exist. 


A very extensive review of subharmonic 
resonance phenomena is given by Ku in his recent 
book4, Actually when one realizes the very 
large Interest shown in subharmonic resonance by 
so many investigators throughout the years, one 
is surprised to find only one reference dis- 
cussing a phenomenon quite similar to the one to 
be taken up for study in the present report. 
Minorsky5 made, in 1955, a study of a nonlinear 
system obeying Van der Pol's differential 
equation. By applying an external harmonic 
excitation to his system he found what he des- 
cribes as an ‘asynchronous action'', It manifests 
itself in the following two ways: 


1. A system subject to limit cycles may 
have its oscillations quenched by the externally 
applied frequency. 


2. A system may be brought into limit 
cycle oscillations by the application of the 
external frequency. 


The word ''asynchronous'!' implies that there 
need not exist any rational ratio between the 
external frequency and the limit-cycle frequency. 
This particular feature brings Minorsky's find- 
ings in close relation to the ones to be 
reported in this paper. There are, however, 
important differences present. Generally the 
above-mentioned asynchronous ''quenching!' and 
"excitation'' actions cannot occur in the same 
system; Minorsky shows that one has to make a 
distinction between the so-called "'soft!' and 
"hard'' systems. Furthermore, he makes no men- 
tion that the system he studied will be neither 
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excited nor quenched by a stationary random 
signal. 


Objective of Present Investigation 


After having demonstrated a simple control 
system that can be excited in an asynchronous 
manner, a general explanation of the phenomenon 
will be given. To understand a phenomenon is to 
have in one's possession a theory with a proven 
capability of prediction. In order to extend 
the ability of prediction the author has not con- 
sidered phase-plane techniques which, although 
used by Minorsky in his above-mentioned work, 
have a very limited range of applicability. 
Instead, an analysis method has been developed 
which permits a study of systems with considera- 
bly higher degrees of complexity. 


As was mentioned earlier, this research was 
actually Initiated as an effort to rid a parti- 
cular system of noise=excited oscillations. For 
this reason the report discusses methods of ob- 
taining most effective suppression of undesired 
oscillations of asynchronous nature. 


The author feels that there might be some 
possibility of utilizing asynchronous action for 
improvements of system response. These thoughts 
are expounded in this paper. 


Theory 


Consider the simple control system depicted 
in Figure 1. Its forward loop incorporates a 


linear (2) and a nonlinear (1) component. The 
linear component has a transfer function: 
K 
G5 (3) = 
n s(1 + s)2 
(1) 
(K = real) 


.The nonlinear component may represent an 
amplifier having a gain-versus~-signal-level 
characteristic as shown in Figure 2. 


Assuming that the system is operated in 
such a way that the error e(t) will not exceed 
€o, i.e., the breakpoint of the gain curve, the 
overall gain will be: 


Ktotal = 


With the above limitation put on the signal 
level it Is clear that the overall system can be 
analyzed as a linear system. In Figure 3 are 
drawn the Nyquist plots for three different 
values of the gain. The system is obviously 
stable for K% 2. The gain is now set equal to 
unity corresponding to a gain margin of 100 
percent and the system response is studied on an 


analog computer. The nonlinearity is simulated 
by means of diodes. 


Figure 4 depicts a series of seven computer 
records. The four simultaneous recordings show 
in this order: 


r(t) = input 
c(t) = response 
e(t) = error signal = r(t) = c(t) 


e,(t) = modified error signal 
(defined in Figure 1) 


Graph 4A indicates the system response for 
a step input and shows clearly the stable 
nature of this system. Without changing the 
amplitude of the step input, the response is 
then studied with a high frequency dither signal 
superimposed upon the input. Graphs 4B, C, D 
and E depict the change in response for increas- 
ing amplitude of this dither signal. Two 
interesting observations are made from these 
graphs: 


1. The dither frequency is so high that 
no trace of it can be found in the system output 
c(t) because of the filtering effect of Gp. 


2. When the dither signal reaches such an 
amplitude that its added effect on r(t) brings 
the nonlinear element beyond its breakpoint, 
(can be seen on the e; recording) the system 
turns more and more oscillatory. 


Should the dither signal be given large 
enough amplitudes, the system will be driven 
into instability as shown in graph 4F. This 
graph indicates actually that instability occurs 
even without a signal input, i.e., it behaves in 
this respect as an unstable linear system where 
oscillations start to build up from a quiescent 
state. Should the dither signal be removed after 
the system has worked itself up to a certain 
amplitude of oscillations, the oscillations 
would either die out or continue to grow. The 
latter would happen if the signal were removed 
so late that the (low frequency) oscillations 
themselves would be large enough to bring the 
error signal sufficiently beyond the breakpoint. 


Graph 4G finally proves what was stated 
earlier, namely that oscillations can be excited 
also by a noise signal. 


The example just demonstrated referred to a 
particular system exhibiting a particular non- 
linearity. In attempting to explain this pheno- 
menon, we should search for a method enabling us 
to preserve a high degree of generality. In 
spite of the difficulty to generalize nonlinear 
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system behaviors, it has been possible to devise 
a technique which permits prediction of the 
system response in the presence of high fre- 
quency dither. The method is an extension of 
the so-called describing function technique. 


The describing function method has been 
proved quite useful in study of nonlinear 
systems. With some modification, this method 
can also be applied to gain an understanding of 
the present phenomena. The keystone of the des- 
cribing function analysis is an assumption, the 
value of which is usually guaranteed by the 
frequency discriminating characteristic of a 
normal type control system. The method fs 
usually capable of predicting not only the 
possible existence of sustained oscillations but 
also the corresponding (approximate) amplitudes 
and frequency. It does not predict the wave- 
shape and is of little value in determining 
transient response. 


The "black box'' N in Figure 5 represents a 
nonlinearity the only restriction on which is 
that it must be non-time varying. The effect 
of this restriction is that any sinusoidal input 
will result in a periodic output the period time 
of which equals the period time of the input. 

We may write the output as a simple Fourier 
series: 


fee) 
y(t) = - By sin (nwt + %,) 
n= 0 


for an input 


(2) 
x(t) =A°* sin wt 


By focusing attention upon only the funda- 
mental frequency component in the output one 
defines a ''describing function'' N, 


° Pal 


By 
A 1 (3) 
which function indicates how the signal changes 
amplitude and phase by passing the ''black box''. 
From then on, linear analysis techniques are 
employed°, One normally has to take into 
account the fact that N changes with both A 
and W, 


After this sketchy review of the describing 
function technique, we return to our actual 
problem. Consider the ''black box'' in Figure 5 
again. Let the input now consist of a sum of 
two separate frequencies: 


(4) 


x(t) = A sin wt + Ag sin Wt 


The output can now be described by a double 


Fourier series/, 
fee) 
Se Bnyn ° Sin C 


n =-00 


y(t) = 
m= 0 
(5) 


(mg + nw)t + nn | 


This statement is not self-evident and 
requires some explanation. A simple heuristic 
proof is given as follows: the nonlinear 
relationship of the "black box!' may be expressed 
as an infinite power series: 


(6) 


This is sometimes referred to as ''curve 
fitting''. Usually very few terms are sufficient 
to accurately describe the nonlinearity. By 
combining Equations (4) and (6) one may obtain 
the output as a function of the input fre- 
quencies: 


fee) 


y(t) = y aj (A sin Wt + Ap Sin Wo ty! (7) 


LO: 


For increasing i-values one obtains (in 
addition to a constant) terms of the following 
types: 

sin Wt 

sin Wot 

sinwt 

2 

sin "Wot 

sin wt * sin Wot 
sindwt 

3 


sin7Wot 


sinzut * sin Wot 
ee. 
sin Wt sin™ Wot 


etcetera 


In other words, the output contains (in 
addition to the d-c component) the following 
frequencies: 
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etcetera 


This is indeed what formula 5 implies. The 
validity of the formula can be verified in the 
case where the nonlinearity is of memory type, 
for instance -- hysteresis. 


The output spectrum is obviously considerably 
more dense than in the case of only one input 
frequency. Computation of the different com- 
ponents is associated with great numerical 
difficulty, even in the case of a simple non- 
linearity. This discouraging fact will not, 
however, prevent us from drawing some general 
conclusions of great importance. 


It should be remembered that in the case 
under study it is generally true that there is 
a considerable difference in magnitude between 
the input frequencies, i.e.: 


1 IR teh e 


Furthermore, it will be stated without proof 
that the coefficients Bm,n in the series 5 are 
of diminishing magnitude with increasing values 
of mandn. (This situation is well-known in 
the case of a simple Fourier expansion.) 
Generally this means that the dominant. fre- 
quency in the lower region of the output 
spectrum is W This very important conclusion 
can be demonstrated for a representative 
numerical case. Assume: 


€ 
u 


10 radians per second 


105 radians per second. 


(= 
i] 


The output spectrum, as obtained from 
formula 5, contains in increasing order of 
magnitude the following angular frequencies 
(d-c component neglected) : 


5 rps (corresponding to m= 1, n = 10) 
10 rps (corresponding to m= 0, n= 1) 
15 rps (corresponding tom=1, n= 11). 


The first and third components, correspond= 
ing to the relative high n values of 10 and ihe 
respectively, may be neglected, and consequently 
the dominant lower component in the Output spec~ 
trum is the one having a frequency equaling the 
lowest input frequency = 10 rps. Corresponding 
amplitude is Bole 

At this point we shall now make the same 
basic assumption as one always makes in order to 
use the describing function technique: due to 
the filtering effect of cascaded components in 
the overall system, all relatively high frequency 
components will be neglected. 


This assumption combined with previous 
findings relating relative magnitudes of spectral 
components gives as final result that the only 
output component of any importance will be Bo ls 


It should be clear at this point that it now 
is possible to define a “describing function"! in 
an analogous fashion as was done in Equation (3): 


B 


(8) 


There is, of course, a fundamental differ- 
ence between the two describing functions defined 
by Equations (3) and (8), respectively. Whereas 
the function defined by Equation (3) depends only 
on W and A, the new describing function wil] 
depend not only on those two variables but also 
upon A, and W,, i.e., the amplitude and frequency 
of the dither signal. For the purpose of iden- 
tification the function defined by formula 8 will 
be referred to as 'modified describing function". 
For most types of nonlinearities, one example 
being the one shown in Figure 2, the dependence 
on the dither amplitude will be pronounced 
indeed, The dither frequency, on the other hand, 
usually has a nonsignificant effect. 


The last two statements are very important. 
In fact, they provide the very clue for a proper 
understanding of not only the asynchronous 
excitation phenomenon, but also the phenomenon 
of signal stabilization. The presence of the 
dither signal controls in effect the gain and 
phase characteristics of the nonlinear element, 
as felt by the low frequency signal. The gain 
and/or the phase of the element may be controlled 
by varying the amplitude but not the frequency of 
the dither signal. This influence of the dither 
amplitude on the low frequency gain and/or phase 
is indeed most interesting. Borrowing an ana- 
logy from the field of chemistry, one might think 
of the dither amplitude as a catalyst, the pres- 
ence of which is necessary for a certain process 
to take place. The independence of the dither 
frequency explains immediately the asynchronous 
characteristics of the phenomenon. The only 
criterion is that this frequency Is high enough. 
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Numerical Evaluation of the 


Modified Describing Function 


So far the discussion has been entirely of 
a qualitative nature. We shall now proceed with 
a presentation of a method aimed at obtaining 
numerical information in any particular case of 
interest. From what was said in a previous sec- 
tion, the effect of dither injection can be 
uniquely predicted on basis of the magnitude and 
phase of the modified describing function N in 
formula 8, knowledge of which presupposes that 


the amplitude By | and the phaselag angle 9%, | 
be known. : : 


The problem of calculating the amplitudes 
B in the double Fourier series 5 should be 
expected to turn out to be considerably more 
difficult than the analogous problem relating to 
a simple Fourier series. The reader is reminded 
that the amplitudes in the latter case can be 
calculated from a simple integral extended over 
the period-time of the function being Fourier 
analyzed. In the present case of a double 
Fourier series, one makes the discouraging dis- 


covery that the amplitude Bm,n must be evaluated 
from double integrals. It would seem at a first 


glance, then, that this would mean an effective 
stop to any further analysis attempt. Fortu- 
nately, this conclusion is not correct. A 
considerable simplification can be achieved by 
once more making use of the assumption that 

Wo” 4. In this case the problem of determining 
B,.1 will be reduced to solving a simple integral. 
This situation will be demonstrated with a simple 
example. Consider again the system depicted in 
Figures | and 2 which was discussed in the first 
part of the previous section. The recordings 
shown in Figure 6 indicate the input to and the 
output from the nonlinear element of this system. 
In this particular case we have 


Wo Gos 


33W 


The output graph clearly indicates a 
periodicity corresponding to the lowest fre- 
quency in the input, i.e., W rps. However, in 
a mathematical sense, no periodicity is present; 
the existence of such a periodicity would require 
that W. were an even multiple of Ww In reality, 
the output contains frequencies below W which 
one realizes immediately from formula 5; one 
can quite easily find integers m and n satis-~ 
fying the equation. 


mW + nw < Ww 


These lower frequencies are indeed present 
in the output wave of Figure 6 but are of such 
a small magnitude that they cannot be detected. 
(The reader Is here reminded of the earlier 
discussion where this situation was explained.) 


It seems intuitively justifiable to 
attempt to calculate Bo,] under the assumption 
that w is the lowest fréquency and thus employ 
the simple formulas of simple Fourier analysis. 
In other words, the output wave of Figure 6 Is 
assumed periodic with a period=-time of 


W 
Lb vis 27 
We obtain 
1 
2 
Boz TY Y(t) * sin wt dt (9) 
fe} 


It can be shown that in general this for- 
mula is quite exact assuming that Wor> Ww, 
Although formula 9 represents the ultimate in 
simplicity that we can expect to achieve, there 
might yet be a big numerical step left to take. 
At this point a digital computer comes in very 
handy. The curves in Figure 7 show in normalized 
units the modified describing function of the 
above-mentioned nonlinearity. The absolute value 
of N is plotted versus dither amplitude (Aj) and 
signal amplitude (A). The phase of N is, of 
course, zero for this type of nonlinearity. A 
few words should be said about the numerical inte- 
gration necessary in this and similar cases. The 
factor y(t) appearing in the integrand will even 
for relatively simple types of nonlinearities 
defy a simple analytical description. This fact 
is readily appreciated after a look at the out- 
put waveform In Figure 6. Each time breakpoint 
of the nonlinearity Is passed, and this occurs 
a large number of times throughout the time 
interval T, a discontinuous shift takes place in 
the derivative of the output. If the integral 
9 were to be solved analytically, one would have 
to first pinpoint all those intervals and then 
perform the overall {ntegration by subdividing 
the integral into a large number of part inte- 
grations. This difficulty is bypassed if a 
digital computer is utilized. In preparing the 
integral for a computer program, it is first 
written as a sum: 


2 
Boe nara y(t) sin wt dts. 


oN 


(10) 


ia T aL 
eae YoCimcn pres in w {it .) 
i=o 


n Is the number of subdivisions of the 
interval 0-T. The choice of n is largely 
dependent on the ratio W/W, because one. must 
make the intervals sufficiently small to 
accurately account for the change in the most 
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rapidly changing term, I.e., be one due to 

sin W.- The program is then worked out {in such 
a way that the computer in sequence computes 
and adds up terms for Increasing values of i. 
In so doing, it should make a logical choice of 
proper value of 


yi -) 


in accordance with: 


If x(t) Sfe,| x(t) 


choose y(t) 


If x(t) > e, choose y(t) = 
=e, t+ ( x(t) - €,) tan % ; 
If x(t)< -e, | choose y(t) = 


=- [eo + ( x(t) -e9) tan a | 


Because of the Inherent capability of the 
computer to make such a choice a numerical 
evaluation of By will be quite simple. The 
curves in Figure’ 5 Peale 35 different com- 
puted values of B - In each case a sub- 
division into 7208’; parts was used. Total 


computing time was 10 minutes on an JBM 709. 


It is interesting to note that the perfor- 
mance predicted from this data was in perfect 
agreement with results obtained on the analog 
computer. 


Effect of Various Types of Nonlinearities 


In the example just studied the presence 
of dither caused an increase of system gain 
with a resulting loss of system stability. 
Other types of nonlinearities might give rise 
to different effects. The saturation type of 
nonlinearity Is quite interesting from several 
points of view. From a historical point of 
view, this type of nonlinearity, depicted in 
Figure 8, was present in a haa where signal 


stabilization first was used The graphs in 
the same figure indicate the modified describ- 


Ing function calculated for this type of an 
element. As might be intuitively felt, the 

gain decreases with both signal level and dither 
amplitude. This is, of course, the explanation 
why Oldenburger was able to detect a stabilizing 
effect in a system incorporating this type of 
nonlinearity. 


In order to further test the validity of 
the theory presented in this paper, a particular 
system was simulated which contained the satu- 
ration type nonlinearity. The system is shown 
in Figure 9 and its response with and without 
dither Injections in Figure 10. Instead of 
obtaining a stabilizing effect as might have 
been expected on basis of Dr. Oldenburger's 
findings, this system obviously goes unstable 


in the presence of dither. Actually, this Is 
exactly what can be predicted on closer study. 
The linear part of the system is namely dyna- 
mically conditionally stable, i.e., it will 
actually be driven into an unstability region 
by decreasing the gain. This fact is easily 
verified by a Nyquist plot. This example brings 
out quite sharply the important point that in 
analyzing the effect of a nonlinearity, one 
definitely must take into proper account the 
linear portions of the system. 


Both of the above-mentioned types of non- 
linearities have real describing functions. It 
is interesting to study the effect of dither In 
a case where the nonlinearity is of memory type, 
i.e., the present output is a function of the 
past history of the input. 


Backlash is a commonly encountered repre- 
sentative for this type of nonlinearity. Dither 
will cause a very interesting change in the 
modified describing function for such an element 
Figure 11 contains a number of plots computed on 
the IBM 709. As in earlier diagrams of the same 
type all symbols are defined in the Figure. The 
graphs reveal that dither will increase the gain, 
In addition, an effect can be noted in the phase, 
The backlash introduces normally a phaselag, the 
magnitude of which increases with the magnitude 
of the backlash (2C in Figure 11). The dither 
will decrease this phaselag; in the case where 
the dither amplitude equals half the magnitude 
of the backlash, the phaselag will be zero. 


This dual effect will quite clearly affect 
the system response. Figure 12 shows a specific 
case. The system is identical with the one in 
Fjgure 1 where the backlash has been introduced 
in the nonlinear "black box'', In this case, 
the dither will obviously improve the stability, 
the reason being that the stabilizing decrease 
of phaselag will dominate over the de-stabi- 
lizing increase of gain. It would be’easy to 
think up systems where the opposite effect 
would be detectable. In Figure 12C Is shown 
the effect on the same system when the dither 
is made up by noise. 


Backlash represents the most complicated 
case for which the modified describing function 
has been computed In this report. It may be of 
interest, therefore, to include the computer 
program (FORTRAN) to give a general idea of 
degree of simplicity (see Figure 13). 


Effect of Random Type Dither 


In the previous theoretical treatment, 
attention has been limited to dither composed 
of only one harmonic component having a fre- 
quency of considerably greater magnitude than 
the cutoff frequency of the control system, It 
was shown that the modulation products could be 
evaluated albeit the numerical effect could 
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become quite strenuous. If several high fre- 
quency components are present, the overall effect 
on the output from the nonlinear element will be 
that more frequencies will be present, i.e., the 
output spectrum will become denser. Due to the 
fact that superposition cannot be employed, the 
problem of numerically evaluating the modulation 
product of interest becomes increasingly diffi- 
cult. For example, should we add just one more 
component to the dither in the previously 
treated, relatively simple nonlinear cases, we 
would find that a numerical computation of the 
modified describing function would be an almost 
impossible task, even using highspeed digital 
computers. 


Random noise, from a frequency viewpoint, 
consists of a continuous spectrum. Even if the 
frequency distribution were known, an attempt to 
find analytically the describing function using 
methods presented in this paper are doomed to 
meet with failure, 


A different approach based on statistical 
considerations was tested but failed to give 
usable results. By ‘'usable'' are meant formulas 
which would render numerical results with a 
reasonable effort. It is the author's opinion 
that this general problem is of such a nature 
that it defies practical computation. 


From a practical point of view it may be 
of more value to report certain semi-empirical 
observations that were made during the progress 
of the research. A large number of analog com- 
puter runs were made in order to test the theory. 
The effect of noise corruption was studied for 
every configuration that was tested. The fol low- 
ing important observations were made: 


1. Whatever the system behavior turned 
out to be in the presence of a single harmonic 
dither signal it remained identical for any type 
of dither waveshape including random noise. 


2. The energy content of the dither sig- 
nal (as measured by its RMS value) rather than 
its waveshape seemed to be determining for over- 
all system behavior. 


It is quite simple to prove by inductive 
reasoning that the first observation has general 
validity. If, for instance, a single harmonic 
dither signal can be shown to have a certain 
effect of the modified describing function, this 
effect will be increased by addition of, say, 
one or more high frequency components. 


It Is obviously of great value to be able 
to predict in a qualitative way what will happen 
in the presence of a noise signal. If we can 
postulate the generality of the second observa~ 
tion, It is also possible to make quantitative 
predictions based on a knowledge of signal 
energy content only. The simplicity of this 


rule is, by the way, appealing from an intuitive 
standpoint. The rule states that one can cal- 
culate the effect on the modified describing 
function of a single harmonic dither signal, 
making full use of the techniques of this paper, 
and then assume that any high frequency signal 
possessing the same energy (RMS value) will 
cause the same change in the describing func- 
tion. Of all possible cases tested on the 
analog computer, this rule gave results that 
were correct within + 5 percent. 


Filtering Considerations 


When the asynchronous action degrades the 
performance of a system and it therefore proves 
desirable to rid the system of it as completely 
as possible, the best method to use is obviously 
the one that removes the cause, i.e., filters 
out the dither from the signal. The obvious 
conclusion that can be drawn from this study is 
that such a filter, in order to be of maximum 
efficiency, should be placed before the non- 
linear element. If the filter were placed after 
this element, it could effectively reduce the 
high frequency portion of the signal, but too 
late to prevent its effect on the low frequency 
gain and phase. The fundamental rule for fil- 
tering would therefore be: reduce the dither 
content to a minimum before it reaches the non- 
linear element. How effective a filter to chose 
can be determined from diagrams of the type 
shown in Figures 7, 8, and 11. 


Asynchronous Action as a Useful Design Tool 


Signal stabilization has been earlier 
mentioned as an example where asynchronous 
action manifests itself in a beneficial way. 
can also be used to obtain an artificial 
in system sensitivity in cases where, for in- 
stance, ''dead-zone'' of backlash result in insen- 
sitive regions for small signal levels. When 
the phenomenon of asynchronous action will be 
better understood, and this paper makes a claim 
to furthering this knowledge considerably, it 
is entirely possible that new uses can be 
envisioned. One feature of asynchronous action 
is the possibility of obtaining a variable gain 
simply by varying the dither amplitude. One 
might expect possible uses of this characteris- 
tic in design of adaptive controllers for 
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instance. [t seems worthwhile also to investi- 
gate compensation networks where one utilizes 
the special features of nonlinear networks as, 
for instance, recently reported by Mishkin and 
Truxal? put combined with the additional 
possibilities of asynchronous action. 


References 


1. MacColl, LeRoy A., Fundamental Theory of 
Servomechanisms, New York, D. Van Nostrand 


Co., Inc., 1945, pp. 78-87. 


Oldenburger, Rufus, ''Sjgnal Stabilization 
of a Control System,'' ASME Transactions, 


Vol. 79, 1957, pp. 1869-1872. 


Oldenburger, Rufus, Lin., C.C., ''Signal 
Stabilization of a Control System,'' AJEE 
Transactions Applications and Industry), 
May 1959, pp. 96-100, 


Ku, YeH., Analysis and Control of Nonlinear 
Systems, New York, The Ronald Press Co., 


1958, pp. 175-190. 


Minorsky, N., ‘On Asynchronous Action,''! 
Journal of the Franklin Institute, Vol. 259, 
195 5 Ppeme209 —=219). 


Truxal, John G., Control System Synthesis, 
New York, McGraw-Hill Book Co., 1955, 


pp. 559-612. 


Bennett, W.R., ''New Results in Calculation 
of Modulations Products,'' Bell System 
Technical Journal, Vol. 12, 1933, pp. 228- 
243. 


Mishkin, E., Truxal, J.G., "Nonlinear Com- 
pensation Networks for Feedback Systems,'' 
IRE National Convention Record, Part 4, 
1957, PPe 37. 


Acknowledgement 


The work reported in this paper was done 
at Hughes Aircraft Company, Culver City, 
California. 


ientl 


1c 


0 3 
Se ds 
w g's o 
gins 
+ Soa om 
o 
are 
BSS 
$733 
FES E. 
oso 8 
nono ft 
gw4H 
aS 
ood 
cee e 
dees 
= 
aaa 
Qo 
=a 
EEE 
: 


(e) indi 


(a) 


1 


ither to be of suff 


F 
e(t), 


in 


e(t), 


ig 
Recordings (f) and (g) prove that oscillat 
order: 


scent state, assuming the d 
t in r(t), 


Response of system shown 
pic 


amplitude 
quie 
dep 


4, 


g 


i 


F 


and e, (t) 


88 


1 


‘omn3Ty oy} ul peyotd 
-op Aj TVsUuTTUOU B IOJ [ADT TeUSTS pue opni[dwe reyIIp 07 
yoodsoed YIM UOTOUN ZUIQIIOSEp PeTJIPOU oY} JO UOTYETACA CYL *L ‘SIA 


a/v 
Oz Si 01 sO. 0 


$ Om uis y+ 4 muis y =(4)x (4)A 


‘g ‘31g Ul UMoYs 
AyireoutTuou eyy Aq Teusis Aouenbe.rj-0M}j B JO UOTYeINPOW °9 “SIT 


1 89 


“AP[IVOUTTUOU JUBTIVAUT-OUNT], “G “STA 


N 
4IN3W313 
YVANIINON 


= x (t) 


x(t)=A sinw t + 
Ag sinwet 


INI 


A,/B 


Fig. 8. The variation of the modified describing function with 
respect to dither amplitude and signal level for a 
saturation-type nonlinearity. 
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Fig. 9. A conditionally stable control system used for simulation. 
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Fig. 11. The variation of the modified describing function (absolute 


value and phase) with respect to dither amplitude and signal 
level for a back-lash-type nonlinearity. 
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— Pagan 


Fig. 12. Response of system containing backlash, (a) without dither, 
(b) with harmonic dither signal, (c) with noise dither. 


DIMENSION C(3), A(6),D(6, 3), FY(6, 3) 


DO8K=1, 3 
DD=K 
C(K)=0.4+DD/ 10. 
DO8J=1,6 

DDD=J 
A(J)=DDD/10.-0.1 
SCM=0. 

SSM =0. 

L=0 

DOYI=1, 1818 

E=I 


X=COS 1F(E*1.5707963/909. )+A(J) *COS1F(E*1.5707963/9. ) 
IF(X-(1.+A(J)-2*C(K))) 1, 2, 2 
2 XX=1.4.A(J)-C(K) 
L=L+1 
IF(I-L)3,3,4 


1 IF((TT-C(K))-X) 4, 4,5 
5 XX=K +C(K) 
GOTO3 
4 XX=TT 
3 YC=(1. /909. )*XX *COSIF(E*1.5707963/909. ) 


YS =(1. /909.) *XX * SINIF(E*1.5707963/909.) 
SCM=SCM+YC 
SSM=SSM+YS 
5) TT = xex 
D(J, K)=SQRT 1F(SCM**2. +SSM**2. ) 
FY(J, K)=(90./1.5707963) *RT ANIF(SSM/SCM) 


8 PRINT 10, C(K), A(J), D(J,K), FY(J, K) 
10 FORMAT (4F 16.8) 
STOP777 


Fig. 13. FORTRAN statement used in computation of curves in 
Fig. 11, 
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ON THE OPTIMUM SYNTHESIS OF SAMPLED DATA MULTIPOLE 
FILTERS WITH RANDOM AND NONRANDOM INPUTS* 


H. C. Hsieh, C. T. Leondes 
Department of Engineering, 
University of California 
Los Angeles 2, California 


Summary 


This paper considers the synthesis of optimum 
sampled data multipole filters with n inputs and m 
outputs. The signal portion of each input is as- 
sumed to consist of a stationary random component 
and a polynomial with unknown coefficients but 
known maximum order. Each signal is corrupted by 
stationary random noise. The filter under invest-— 
igation is linear, time-invariant, and has finite 
memory. Each input to the filter consists of a 
sequence of impulses with a constant period T. 
Each impulse is assumed to have an area equal to 
the value of the signal plus noise at the sampling 
instant. 


The synthesis procedure to be developed is to 
specify the weighting functions of the filter such 
that the system error, which is defined as the dif- 
ference between the actual and ideal outputs, has 
zero ensemble mean and the system ensemble mean 
square error is minimum. The weighting functions 
thus obtained will have, in general, abrupt jumps 
at the sampling instants but they are continuous 
within the sampling intervals. 


The synthesis procedure is extended to the 
case shown in Appendix A where each of the nonran- 
dom signals can be expressed as an arbitrary linear 
combination of a set of known time functions. Fur- 
ther generalization is possible to the synthesis 
of time-varying filter with sampled nonstationary 
random inputs as given in Appendix B. 


Introduction 


A Survey of Two-pole Linear Filtering Theory 


A two-pole system is defined as a device which 
has only one input terminal and one output termiml. 
During the past few years, considerable attention 
is given to the synthesis of filters whose signal 
contains both randop and nonrandom components. 
Zadeh and Ragazzini~ have solved the problem of 
synthesizing optimum, continuous, linear, time in- 
variant filters, having finite memory time, for use 
with signals consisting of a stationary random pro- 
cess and a nonrandom polynomial with unknown coef- 
ficients but known maximum degree n. They assume 
the signal to be corrupted by an additive station- 
ary random noise. The,Zadeh-Ragazzini filter can 
be reduced to a Wiener~ filter if the nonrandom 
part of the input is taken to be zero and the mem- 
ory time is infinite. It has been shown that the 
optimum weighting function will involve impulse 
functions of various orders. 

*This research was supported by the United States 


Air Force under Contract No. AF 9(638)-438 moni- 
tored by the Air Force Office Of Scientific Re- 


search of the Air Research and Development Command. 
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With the development of sampled data systems, 
statistical design techniques have been adapted to 
the case where the input is pogertce only at dis- 
crete instants of time. Franklin- considered the 
optimum synthesis of a sampled data filter and on 
trol system with input being the sum of a station- 
ary random message and a stationary random noise. 
His problem is exactly the sampled data counter- 
part of Wiener's theory. For the optimm filter 
the transfer function is obtained by using the us- 
ual spectral shaping method. It will contain, in 
most cases, the product of two terms, one of which 
is a rational function of S and the other is a ra- 
tional function of e**. The control system prob- 
lem has also been investigated by Chang+. However, 
his synthesis procedure is carried out mostly in 
the Z-domain as contrasted to Franklin's time- 
domain approach. 


The extension of the optimum filtering theory 
of Zadeh and Ragazzini_to the sampled data case 
has been given by Lees’. His optimum linear fil- 
ter will convert its input sequence of impulse 
functions into a smoothed output subject to the 
following conditions: the weighting function of 
the filter is not zero over a finite range only; 
in absence of random components of the input, the 
interpolation or extrapolation is error-free; in 
the presence of any random signal, the mean-square 
error of the output is minimized. It has been 
shown that the weighting function of this optimum 
filter will be piecewise continuous within the 
sampling intervals. Later on, Bergen- also con- 
sidered the optimum linear filtering of a polynom- 
ial message plus stationary random noise. His in- 
vestigation is very similar to that of Lees but 
his approach is more conventional. 


Recent Development of Multipole Systems Theory 


As stated in the first section, the t wo-pole 
filter problem has been intensively studied by 
previous investigators. Not until recently has 
the synthesis of mltipole filters or control. S43 
tems received due consideration. The authors’? 
have developed the solution of optimum weighting 
functions in the Wiener sense for the multipole 
filter or control system where the input to each 
terminal consists of stationary random signal and 
noise. A set of integral equations is obtained by 
minimizing the mean-square errors of the system 
outputs. Through transformation, these integral 
equations can be converted into algebraic equa- 
tions. For filter or free configuration control 
system, the optimum system transfer function can 
be found by solving these algebraic equations and 


using the method of undetermined coefficients. In 
the semi-free configuration control system, the 


optimum transfer functions of compensation can be 
obtained in a similar manner such that the outputs 
of the fixed plant are more desirable. 


With the increasing occurrence in practical 
problems of inputs available only at discrete in- 
stants of time, it is highly necessary to develop 
a theory for the optimum synthesis of sampled data 
multipole filter or free configuration control sys- 
tem. This paper presents the theoretical design 
equations for this important problen. 


The Linear Filtering For Sampled 
Data Multipole Systems 


Statement of the Filtering Problem 


The filter under investigation is a sampled 
data multipole system with n inputs and m outputs. 
In order to facilitate the formulation of the prob- 
lem, sampling switches are fictitiously inserted at 
the input terminals. The sampling process is then 
obtained by closing these switches briefly every T 
seconds such that the inputs are observed at these 
discrete instants of tine. 


The input I,(t) to each terminal before samp- 
ling is continuous time function and will consist 
of signal and noise with the following properties: 


(1) The signal S;(t) has two components. One 
of them is a stationary random process M,(t) and 
the other one is a nonrandom function P,(t) which 
can be represented by a polynomial of known maxi- 
mum degree m, but with arbitrarily unknown coeffi- 
cients. 


(2) The noise N,(t) is also a stationary ran 
dom process. 


(3) It is assumed that the statistical char- 
acteristics of the random signals and noises are 
completely known. 


With the above statements, the total signal to 
each terminal of the filter can be expressed by 


S(t) = M(t) + Py(t) (1) 


where 
; a ti (2) 


P;,(t) = 
Tr=0 


SMeyy te ea IS CAS Ba sl 
The input I,(t) is then 
Iy(t) = M(t) + P(t) + H(t) (3) 


It should be noted that the orders of the polynom- 
ial signals P,(t) can be different from terminal to 
terminal in order to give the best representation 
for the nonrandom inputs. 


Because of the extremely short closing dura- 
tion for the samples in comparison with the samp- 
ling interval T under actual physical situation 
and also because of the resultant mathematical 
simplifications, the sampling operation is taken as 


an inpulse modulation of the continuous function. 
Thus the sampled input is 


T(t) = I(t) Saal 5 (t-i?) (4) 


=— 09 


where d (t) is an impulse function. The system 
block diagram is shown in Fig. 1. Here the W5x'S 
are the weighting functions between the various 
input and output terminals. 


The synthesis procedure to be developed is to 
specify all the system weighting functions satis— 
fying the following conditions: 


(1) Each system error, which is defined as 
the difference between the actual and ideal out-— 
puts, has zero ensemble mean; i.e. 


E x(t) = 0 (5) 
fore je=) jae, «own 
(2) Each ensemble mean square error of the 


system is a minimun with respect to the weighting 
functions Wj,'s which satisfy the condition l. 


2 ; t 
5 (t) = min. w.r-t. Wy,'S (6) 


LOT? Vjesn Ly, 2; eoe M 


(3) This miltipole filter is physically re- 
alizable, linear, time invariant, and with finite 
memory. The last statement implies that the 
weighting functions W4,'s should vanish for all 
time instants greater than a certain specified nm- 
ber. For simplicity, this number is chosen to be 
a multiple of the sampling interval T. Thus we 
have 


W5x(t) =a for t <0 (7) 


and t > (N41)? 


where N is an integer. 


It is evident that the synthesis approach is 
first to obtain the constraint equations imposed by 
condition (1) and (3), and then to find the optim- 
um weighting functions which will minimize the sys- 
tem error subjected to these constraint equations. 
The filter is consequently considered to be the 
best under these design criteria. 


Generation of System Error 


The filtering problem to be discussed in this 
report is in a very broad sense where the ideal 
outputs can be any linear operation on the signal 
portions of the inputs. Let the ideal weighting 
functions for the system be Dj_(t)- Then the de- 
sired outputs are 


Dy 7) S(t -T) aT (8) 


— oc 


n 
C(t) = ae 
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The lower integration limit of "-oce" emphasizes 
the fact that ideal weighting functions can be phy- 
Sically unrealizable. 


The actual filter outputs can be expressed as 


to 


if W567) Tt -T ) aT 


° 


n 


(9)* 


where Ws,.(7) = 0 forT < 0 and T >t,. By de- 
finitioh, the system errors are 


E(t) 


iT] 


C4(t) - Cqs(t) 
t 


fe} 
n 
Ett ot.) a7. 
os [xv i ) 


° 


n co 
xe [>,7) S(t -T) aT (10) 
k=) 


—_- oo 


The formulation of this problem will be best 
understood by referring to Fig. 2. For simplicity, 
only a single line flow diagram is drawn. It is 
apparent that the system outputs are compared with 
ideal outputs at all time so that the filter will 
be a continuous one. 


Derivation of Constraint Equations 


The constraint equations shall be obtained by 
applying the condition that the mean system error 
is zero. Combining Eqs. (1), (3), (lL) and (10) 
gives the system error as 

t 
° 


n O°, 
Gt) = Ze [Ru T ET S (t-T-it)ar 
° 


n oo 


ts 
n 
= oa i My T[ Wlt=-T) + Ht 7] 
° 


x 22 § (t-T-iT)dT 


*the case where t, differs for each input when a 
particular output is under consideration, is a re- 
latively straightforward extension of the work 
presented here. 
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(11) 


If it is assumed that the ensemble averages of 
the random signals and noises are zero 


HAT) = HIE = 0 


it follows that 


(12) 


1=-a@ 


n (2 aad 
€;(t) = | | Wil TP (t-T) pa (t- T-iT aT 
fe] 


is f Dy T )P (t= T)AT (13) 
a) 
With x169) = 0, the following equation is gen- 
erally obtained. 
te oo 
| Wye T Py (t= T a S (t-T-iT)aT 
° 
= J Dsy( T )Py(t-T aT (1h) 
= 40 


for ke= se 50666 D: 


It is evident that, with zero means for the 
random components of the inputs, the first condi- 
tion imposed on the weighting functions can be re- 
phrased as: In the absence of random components 
of the inputs, the actual filter outputs are iden- 
tical to the ideal filter outputs. 


The function P,(t-7) is now expanded into 


Taylor's series about Py(t). Since P,(t) is a 
polynomial of m, degree, this expansion will give 


Ti r 
Py(t-T) = os a Perley TF (15) 


By substituting Eq. (15) into Eq. (1h), there 


results 

to 

me yer r os 

pia {ou P, (+) [7 WT) 2 S(t-T-in)ar 

r=0 . i=-60 
) 

My, - = (r) ri 
= ee P,. (t) fr Ds,(T aT (16) 


If this equality is true for all values of r, the 
following condition must exist. 


us Pa. 
| Ty (T) 2 § (t- T-it) aT 
1=—60 
fe) 
a r 
2 i a Dy TAT (17) 
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The left side of Eq. (17) needs further inves- 
tigation. Since W;,(T) is also zero for T >t,, 
the range for the nde iis restricted. It can 
easily be verified that 


t—-t 
co) 


(18) 


Define 


iz=— oo 


t oo 
(e) 
Ur = J y W5y(T ) De S (t= T-iT)aT (19) 
fe} 


Then, by equation (18), it can be expressed in sum 
mation form as 


“ 

iS 

2 int 

se (t-a9)” Wy (t-i7) 


Ue = (20) 
i= 
Wy 
Let t be replaced by 
tos hatha) (21) 
where h is an integer andOQG A =< T. Then Eq. 


(20) can be rewritten, with the further substitu- 
tion of t, = NT +A, as 


aly 


ots > (hT-aT +A)" Wy (hT-iT +A) (22) 
i=h-N 


n+ & 


If now h-i is replaced by q which must be an in- 
teger, it can easily be shown that 


N 
Up = ys (q? +A)" Wy, (aT + A) (23) 


Q=o 


The range for the index q is quite consistent with 
the condition that W5u(t) = Ol fort —Orand 

t > (N41)T as it should be. It is evident that 
the weighting function W.,(t) is defined piecewise 
during each sampling interval. 


The right side of Eq. (17) may be denoted by 


ve if a" Dal TAT (2ls) 


=<sd0 
where Vig is a constant. 


Then Eq. (17) can be stated more conveniently by 


te =i (25) 


igor Wis SS Bh By ago tl 


i EO). ike coe my, 


n 

These are 28 m, + n constraint equations which 
k=) 

W5x(qT +24) must satisfy for all values of A . 


By inspection of Eqs. (23) and (24) it may be 
noted that, in order to obtain error-free outputs 
in the absence of random components of the inputs 
and to minimize the mean square errors in the pre- 
sence of them, the number N for the filter must be 
greater than the maximum order of all the polynom- 
ial inputs. In other words, the filter must be 
able to remember more than m, +1 previous sample 
values. 


Optimum Synthesis of Sampled 
Data Multipole Filter 


Minimization of the Mean Square Error 


Suppose now that the constraints given by Eq. 
(25) are satisfied, the expression for error will 
be reduced to 


t 


ie) 
n 
G(t) = me J WT) | alt Tame 7] 


te) 


Ke GS (t-Taar 


j=-<0 


n iv.) 
= 2 | Dal Tilt Tat (26) 
re gk\ 7 Me 
=— 29 
Let t be replaced by hT + €,as before. Then, by 
using the same manipulation as in Eqs. (22) and 
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(23), it 
E(t) 


is easily shown that 


iT] 


€5 (ar +2) 


n [o2) 
= Mes WT) | Myce T)04(-7) | 
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“ 2 i Dal T (tT) T 


The mean square error is then given by 


(27) 
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(28) 


where Br and @,,, are the summation of various 
autocorrelation functions and cross-correlation 
functions for random signals and noises as defined 


by 


BT) = By gu, Fu, nT Fy uD Oy nM 
(29) 
O4(T) = Ay vu,( + by u,(T (30) 
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The above correlation functions are the ensemble 
averages of the random variables. If the ergodic 
property is applicable to these random processes, 
these ensemble averages can be equated to the av- 
erages with respect to time performed on a single 
representative function of the ensemble. 


The next step in the development is to mini- 
mize the mean square error with respect to the 
weighting functions subjected to the constraints 
in Eq. (25). By introducing the Lagrangian multi- 
plier Cie » this problem is equivalent to minimiz—- 
ing the peeneicc 


mn Mr 
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- 09 
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— 00 
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(31) 


The technique! 2 employed in minimizing Eq. 
(31) is essentially to complete the square for the 
first two summation terms. Let the arbitrary func- 
tions Us,(T)'s be defined implicitly by the fol- 
lowing dane tans’ 


n oo N 
= U5 Tr )Brjcl To- TZ AT --A)a Ty 


(e) 


Mi r 


n oo 
= 2E i Dail TO x0 46 T o> Ty)d aye, Ak'rT2 


(32) 


By substituting Eq. (32) into Eq. (31) and adding 
and substracting a term 
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(33) 


It is quite evident that the first three sun- 
mation terms can be expressed as 


ape 
oe i [ mT) T)]f M(t T Ait T)] 
fe} 


—_—_--—_--—---— 5 


N 
x2. §( T-of-A)aT (3h) 
q=o0 


which shows the fact that this term is always pos- 
itive. Thus, to minimize Eq. (33), it is obvious 
that the following choice can be made 


U5x(T ) = Wail T ) (35) 
Under this condition, Eq. (32) becomes 
n oa N 
Wa (Ty Br (To-F ) T-4)d 
z. | gc Prof meses: aT, 
ro) 
n Mt 
= r 
= i Pil TIO 14¢6 T 9> TAT} ye Meir T 2 
forukte xtc yecec 
To = AT +A where -=0, 1, ...¥ 
(36) 


To show that Eq. (36) gives not only the ne- 
cessary condition for optimum system but also the 
sufficient one, let the optimum weighting functiors 
W.,(T ) be replaced by another arbitrary weighting 
tinction of the form W5x(T) - Y5x(T )- Then Eq. 
(31) becomes 


n_ 00 N 
= 25 + 2 Bite Te) ZS {1-4 )dT, 
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[Efe 


n °o N 
a aE { Maid Ty Tor TZ ST 4) art. 
° 


Myr 


BACT Deyn (s- 1147 - me AwirTs 


eo 


n 
rT: tn | sin || 
+ | Y 5x6 fa )+ N,(t-= T) 
2 
N 
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(37) 


Thus, if Eq. (36) is satisfied, z, reduces to 


1 = ma 
Z.=2;+ 2 | tx] wT nce) 
ce) 


LE § (T-qt- a)aT 


g=o 


(38) 


Hence, it is seen that the value of Zs; is not de- 
creased by any perturbation of W., if*Eq. (36) is 
fulfilled. The sufficient condition is thus es- 
tablished. 


Solution for the Optimum Weight Functions 


In order to solve the optimum weighting func- 
tions, Eq. (36) is rewritten for investigation. 
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Eq. (39) can be rewritten in summation form as 
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(41) 
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These are n x (N +1) simultaneous equations 
which can be used to solve for MOR +A) in 


terms of A kip? By using the os m, +n constraints 
given by Eq. (25), these Lagrangian multipliers 

will be determined. Thus, combination of Eqs 
(2b) and (1) gives a unique solution for the op- 
timum weighting functions of the filter. 


For convenience, let Eqs. (1) and (25) be 
rewritten as 
my, 


os = Be ry( LI-a0)W yy (GTA) oe eects 
k=1 q=o 


= By £1 +d) (424) 


Imfoy Mell KS ae Lo Son el 


BY Sacre 


and 
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fore kes, Ps eoe Nh 
r=0, ibs evo Mm, 


In order to solve Eq. (2) by matrix algebra, 
it will be desirable to assume that the polynomial 
inputs have the same order ap For any polynomial 
input which has order less t an m,, its correspon 
ing Lagrangian multipliers Ayu; "Pin subscript 
r > mr will be zero. In addition, Eq. (2B) will 
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be identically zero for these subscripts. Thus 
the total equations to be considered now are 


ax [+2 +( +0]. 


Let the following submatrices be defined: 


jk! 


Bsyr(Q) 
Bar (T+ J) 


B yy (NT A) 


W5x( 4) 
W5y(T+ dd) 


= (N+1) column vector 


= (m,+1) column vector 
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= (N+1) column vector 


= (m,+1) column vector 
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Ao LE eA 
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— 5 : = (m4 
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Ale (ts A)eeo et (nT+4) 
J ate = : f = (m+ 
ecake : 3 5 x(Ne) 
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(U9) 


AS pire noe 2 (wr+a) P 


It should be noted that, with the actual order 
of any polynomial input less than the maximum order 
» the appropriate elements in submatrices G skh 

Kt» Ips and V5, will be zero. ms 
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Then Eq. (2) can be rewritten in normalized matrix 


form as 91D eeeeeee O 
Hyy Hy eeee- Hy! : , , QO do ceeeeee 0 


Hoy Hoo ecoeoe | 


Ie 
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ci : (53) 


Jn 
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2 | 
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Eq. (50) will be expressed in a more compact _jn 
form by further defining 
E; 
1 
Hii the eeeeevee Hin UJ 
H bie 
Ho, Hoo eeeeeee en é 
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Ty ie) esceeese @) 32 
O Ty eevee 0 Vj =| + | = column vector (58) 
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: = (59) 
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Eq. (59) may be solved to yield 


F,| (2 8 | [3 
= (60) 
ay] |Z 21] 
where 
Riese Bsn Wiest (siinat)id/ te) (61) 
8s Hi r(se? 1)” (62) 
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othe fie (64) 
Hence the weighting functions vector is given by 
By Roh Ey +35 (65) 


imore @) Stal. (5 eoe M 


It is evident that the submatrices R and S are the 
same for all output terminals. However, the E 

and V4 submatrices will be a function of the — 
ticuter output terminal concerned. 


Although Eq. (65) gives a concise expression 
for the optimum weighting functions matrix, yet 
the evaluation of R and S as given by Eqs. (61) 
and (62) may actually be very tedious. It can be 
seen that these functions have, in general, abrupt 
jumps at the sampling instants but they are con- 
tinuous within the sampling intervals. If A is 
equal to zero, a discrete filter is obtained. 


Mean Square Error of Optimum Filter 


The mean square error of a particular output 
terminal for the filter is given by Eq. (28). It 
can be rewritten as 
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Since the 


where B.,,(7) is defined by Eq. (0). 
random dEenals and noises are stationary and the 
sampled inputs are non-stationary, the mean square 
error is a function ofA . An average of the mean 
square error over O< A <_ T will yield the av- 
erage over all time. 


The minimum mean square error will be obtained 
if the solution for optimum weighting functions as 
given by Eq. (65) is substituted in Eq. (66). By 
using the matrix notation, it can be conveniently 
expressed as 


ii] 

ty 
ifs} 

ie) 


2 
AR Nae, 


x(R BE} +S V3) +45 (67) 


Here the prime denotes the transpose and Aj is de- 
fined by 


Voor oz 


Ey Dy Te8T 


ao 


«( D5x( TF a! To- TdT, 


= 


(68) 


An Illustrative Example 


As a simple example, let a filter with two in 
put-terminals and one output—terminal be consider- 
ed. The signals to these two channels are assumed 
to be the same. However, they are corrupted by 
different noises. The polynomial inpus are taken 
to be of first order. The random signals and noi- 
ses are not correlated. Furthermore, the noises 
are uncorrelated through the sampling interval T. 
The function of the filter is to extract the sig- 
nal from these two different channels. The system 
block diagram is shown in Fig. 3. 


It is evident that, in this problem, 

M(t) = M(t) = M(t) (69) 

Py(t) = Po(t) = P(t) = A, + Ayt (70) 
Let the autocorrelation function of the random 
signal M(t) be 

ee cake 

A(T) =e (71) 

Under the assumption of uncorrelated noise samples 


through the interval T, the two autocorrelation 
functions for noises are given by 


A = * 
Ay w,{ T-qT) = 0.56 oT (72) 
(At-aT) = 0.2 (73) 
AN, -q ry ee 


where 5 L is Kronecker delta which is either one 


or zero according as the subscripts 4 and q are or 
are not the same. In addition, 


Byy (T) = By (T) = 0 (74) 
1 2 
The sampling interval is taken to be one second and 
the memory of the filter, three seconds, (N = 2). 


The ideal weighting functions of the filter 
are 


D4(T) = (T) = S(T) (75) 


The overall correlation functions defined by Eqs. 
(29) and (30) are 


,(T) = Bo3(T) = Dry T) 


BT) = thy T) + dyin, (T) 


203 


Q(T) = Oye(T) = On(T) = O,,(7) = A(T) 
Let the following submatrices be formed: 
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Even in this simple example, it is difficult to 
calculate the submatrices R and S interms of A 
as given by Eqs. (61) and (62). 
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gE? I 


iT] 
jaa) 

L} 
te2} 
rs 
=~ 
Q 
x 


However, by using digital computer, the numerical 
values for the weighting functions can be easily 
obtained for every assumed value of Q\. 


Conclusion 


A synthesis procedure has been developed for 
obtaining the optimum sampled data multipole sys- 
tem with random and nonrandom inputs. With the 
assumption that the nonrandom signal components 
can be represented by polynomials, the solution 
for these optimum weighting functions becomes a 
simple routine matter. Theoretically, an analytic 
expression can be obtained for each weighting funo- 
tion in terms of A over each sampling interval T. 
However, due to the inherent complexity of multi- 
pole system, the matrix operation for obtaining 
the solution becomes very tedious even in the sinp- 
lest case if A is carried all the way through as 
an independent variable. Thus digital computer is 
advised to be used so that the numerical answers 
can easily be found for every assumed value of A . 
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Appendix A 
Filtering With Generalized Nonrandom Inputs 


In the previous discussion, assumption is 
made that the nonrandom signals P;,(t) can be ex- 
pressed by polynomials with known maximum degrees 
m, but unknown coefficients A,,. It is seen that 
the function P;,(t-7 ) soandadeints Taylor's ser- 
ies about P,(t) as given by Eq. (15) is a poly- 
nomial in 7 with maximum degree of m, and coeffi- 
cients as functions of t. Hence Vy, obtained from 
Eq. (2h) are constants. The continuous filter 
thus obtained is time invariant. 


Extension of the solution can be made to the 
case where the nonrandom signals are not restrict- 
ed to polynomials of degree m, but can be expres- 
sed as an arbitrary linear combination of a set of 
known time functions. Thus, the nonrandom signals 
are given by 


mn, 
P(t) = Sa Alp Pry (t) 


r=o0 


(A-1) 


where p,,(t) are known but A), are unknown. For a 
considerably large class of functions py,(t), they 
have the property that p,,(t-7) can be represen- 
ted in the form 


Mir ys + 
Pyr(t-T) = reset) bt) (4-2) 


Under this condition, the constraint equations are 
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The simultaneous equations obtained in mini- 
mizing the mean square error are 
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It is evident that, except slight modification of 
the matrices I and J, the techniques developed pre- 
viously are equally applicable to this generalized 
nonrandom inputs problem. 


Appendix B 


Time-V Filters With Sampled 
Nonstationary Random Inputs 


The optimization theory thus far developed is 
based upon the conditions that the random inputs 
are stationary and the ideal weight functions are 
time-invariant. During recent years, the analysis 
and synthesis of time-variant systems have receiv- 
ed considerable attention such as those systems in 
the guided-missile problem. Thus the desired out- 
puts of the mltipole filter can be time-varying 
linear operation upon the signal components of the 
inputs. In many practical problems, the assump- 
tion that the random inputs are stationary is not 
valid. The nonstationary random inputs imply that 
their statistical properties varying with time. 
Therefore, the synthesis of time varying filters 
with sampled nonstationary random inputs need to 
be considered. 


For this more general problem, the nonrandom 
inputs will be assumed of the form expressed by 
Eq. (A-1). 


wD 
P,(t) Ss Aur Pr(t) 


r=o 


The desired outputs are then 


n co 
Cg5(t) = 2 fetes | meCt- TP tT) dT 
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(B-1) 


Here a time variable weighting function h(t,7T) is 
defined as the response of a system at time t to a 
unit impulse applied at time t-T, i.e., YT time- 
units earlier. 


The system errors are expressed by 
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For €.(t) = 0, the constraint equations are, with 
t=hT+Q, 
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If these constraints on W(t, 7) are satisfied 
the mean square errors are 
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By using the same minimization technique to 
the system ensemble mean square error, the simul- 
taneous equations which the optimum weighting 
functions must be satisfied subjected to the con- 
straints can be shown as 
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6  (ht- £T ,hT-q? )W., (hT+ A ,qT+ 2) 
es Ps k'k jk 


as L(t)=M+PEN, 
= = : = pecs ie Va ee i 
= By (at+ 4, £044 ) a antsPacrn(Bt- 21) 
(B-5) 
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L,(t)=M+P#Nz 


In order to use matrix algebra to solve Eqs. (B-3) 
and (B-5), a set of submatrices H, I, J, Fy, Gs, 

E; and Vs can also be defined. The najortitrer- les Filth 

Bice from the previous treatment is that the ma- Ale 3 2+ foles Filter 
trices now are a function of t (or h) as well as 
4. Except for the simplest problem, use of dig- 
ital computer is advisable in obtaining the solu- 
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FIG. 2 Error Generition Diagram 
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Fig. 4(a) Weighting function for the example. 


Fig. 4(b) Weighting function for the example. 


207 


{e) 
0 0.1 02 03 0.4 Oo5 06 0.7 O8 O9 1.0- 


Fig. 5. Minimum mean square error. 
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DECOUPLING TECHNIQUES IN MULTILOOP CONTROL SYSTEMS 


R. H. Loomis 
Westinghouse Electric Corporation 
Baltimore, Maryland 


Summa ry 


The stability and accuracy of a control 
System consisting of two interwoven loops are 
affected by interaction among all transfer 
functions and all inputs. This is undesirable, 
since the design of each transfer function must be 
compromised to account for the interaction 
effects. 


In some cases, the two loops can be 
decoupled by the addition of an appropriate feed- 
back path. This paper describes the application 
of this technique to a radar angle tracking loop. 


The technique described has the following 
advantages: 


i The required matching takes place 
within the loop so that the effect of any mismatch 
is attenuated by the open loop gain. 


(be For the example used, the output 
angular rate accuracy becomes independent of 
aircraft motion. 


33 Subject to the rather minor 
restriction implied by item 1, the two loops 
become independent as far as stability is 
concerned. 


Decoupling, which is described in this 
paper by illustrating its application toa 
particular system, is applicable to any complex 
control system. However, tothe author's 
knowledge, no formal procedure for general 
application is available. 


Introduction 


Multiloop Control Systems 


This paper outlines a technique whereby the 
several interwoven loops of a complex control 
system can be made to act, dynamically, as if 
they were independent and operating in cascade. 
This technique can be useful in improving system 
stability in cases where the transfer functions 
are fixed by other considerations. It is 
especially useful in cases where the controlled 
variable is not of principal concern(e.g., when 
output rate accuracy is of principal interest) and 
has been successfully applied to a system having 
this particular characteristic. 
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The control loop shown in block diagram 
form as figure 1 is typical of the class of loops 
to which decoupling can be usefully applied. 


The outer or control loop serves to track 
the reference variable, R, while the inner or 
isolation loop acts to prevent system response to 
the disturbance function, U. However, since 
they are interwoven, both the control and 
isolation characteristics will be affected by both 
transfer functions. For example, the outer feed- 
back path (through G,) can act to improve the 
system isolation capability and decrease the 
phase margin. 


Decoupling 


The destabilizing effect of the outer loop can 
be counteracted by appropriate modification of the 
inner loop transfer function. However, if this 
approach is used, the added isolation brought 
about by the outer loop is lost as the design of 
each loop must be compromised to account for the 
effect of the other loop. 


An alternate approach is to change the 
tracking configuration in order to prevent inter- 
action between the loops. This can be done quite 
easily, although at a sacrifice in control 
"tightness.'' (That is, the output variable (C) 
will not follow the input variable (R) as well.) 


In the example discussed later, the effect 
of decoupling is to make the system shown in 
figure 2 dynamically equivalent to that system 
shown in figure 3. 


The use of a decoupling loop introduces 
problems of drift and component matching. The 
component matching problems are not significant, 
since the matching is done within the loop and, 
therefore, the effects of any mismatch are 
attenuated by loop gain. The drift problem, 
however, requires careful attention to component 
selection. 


Conditions 


Decoupling is not always a desirable 
approach to a stability problem, since the 
increased stability is obtained at the expense of 
loop ''tightness.'' It should be used only after a 


study of accuracy requirements has shown that 
the resultant loss of accuracy is acceptable. 


The following conditions are examples of 
the factors which must be present to justify 
consideration of decoupling: 


1 A system of the general type shown 
in figure 2 is to be used. More specifically, a 
complex control system is to be used. 


ae The transfer functions cannot be 
designed to meet all requirements and still retain 
adequate overall phase margin. 


3. Some degradation in output accuracy 
is acceptable. 


A situation which meets these conditions 
is discussed under ''Practical Use of Decoupling." 
The application of the technique is as described 
under ''Application of Decoupling. "' 


In addition to the above usage criteria, the 
application of decoupling places requirements on 
all components in the loops involved. Although 
the principal purpose of this paper is to present 
an analytical description of decoupling 
techniques, the successful application of these 
techniques is so heavily dependent on the use of 
appropriate components that a brief discussion 
of the component selection problem is included 
under ''Component Matching" and ''Component 
Drift. 


Application of Decoupling 


System Configuration 


The control system shown as figure 2 is 
typical of the class to which decoupling can be 
applied effectively, and it can be used to 
illustrate the decoupling technique. For clarity, 
a specific open loop transfer function will be 
assumed. Figure 4 is an open loop gain 
frequency plot of the assumed function. 


Original Loop Stability 


In the assumed control system (figure 2), 
the outer loop can be regarded as a feedback path 
around the isolation loop for stability analysis. 
This has been done for a particular group of 
parameters, and the result is shown in figure 4. 


Examination of the figure shows that the 
control loop feedback has two effects (other than 
its primary tracking function). First, it 
increases the isolation performance of the loop, 
and second, it decreases the overall phase 
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margin. For the assumed parameters, the 
decrease in phase margin is such that the overall 
loop phase margin is reduced to an unacceptable 
value. 


The magnitude of this reduction is calcu- 
lated as follows: 
ipl SS fe (1) 


Adding the control loop, the total feedback 
becomes: 


H'=s+G, (2) 


Taking the ratio between the two feedbacks 
and using assumed parameter values give: 
s+ G, G 


GH' 1 20 
ae == (1+8/G,) == (148/20) (3) 


The outer loop introduces an increase in 
slope below a frequency of 20 radians per second. 
This is illustrated in figure 4 where the solid 
lines show the stabilization loop alone and the 
dashed lines show the overall system. The 
destabilizing effect is apparent from the figure. 
For these parameters, the loss in phase margin 
is: 


Fane 6) 


90° - tan 30 


23S (4) 

This loss of 33 degrees reduces the overall 
phase margin to approximately 10 degrees, which 
is not a satisfactory value. 


Original Loop Performance 


The performance of the original system can 
be obtained by direct solution using the block 
diagram shown as figure 2. Thus 


s/G)G 


Ce : R+ U (5) 
1+8/G,+s/G,G, 1+8/G,+8/G,G, 


The added isolation from G, is apparent 
from the equation. Figure 10 shows the response 
of this loop to a severe transient in U. From this 
figure, it is apparent that excellent positional 
isolation, as well as low phase margin, is 
afforded. 


An obvious method of improving the phase 
margin is to add a compensating network to G>. 
However, this approach would also eliminate the 
added isolation brought about by the outer loop 
and, in some cases, redesign of the loop would be 
severely limited. (For example, the quadratic 
indicated in figure 4 might be due to a mechanical 
resonance beyond the designers control.) This 
being the case and taking into consideration the 
conditions of a previous section, decoupling may 


offer a more promising approach. 


De coupling Loop 


The control system of figure 2 can be 
decoupled by adding a feedback path from the 
output of the integrator to the input toG_,. In 
general, this would be quite simple to mechanize. 
Figure 5 is a block diagram of the decoupled 
configuration. 


Decoupled Loop Stability 


The effect of decoupling on stability can 
best be examined by writing the control equations 


for the loop. This gives: 
G,/s G, 1 
iE = 
G78) ‘ne! ees TG,) (6) 


Equation 6 can now be used to derive the 
dynamic equivalent of the decoupled configuration. 
Figure 6, a block diagram of this equivalent 
system, shows clearly that the effect on stability 
of loop interaction has been completely eliminated. 
The addition of the decoupling loop makes the 
resultant system act as if it were comprised of 
two independent loops operating in cascade. This 
has obvious advantages in design. Since the loops 
are independent as far as stability is concerned, 
a considerably wider choice of parameters is 
available to the designer. However, for the 
purpose of this paper, the transfer functions will 
be considered fixed and the effect of decoupling 
the original loops evaluated. 


Decoupled Loop Performance 


Equation 6 indicates the effect of decoupling 
on the accuracy and isolation capability of the 
system. If this equation is rewritten, it can be 
put in a form more suitable for comparison with 
equation 5 which is the equivalent relation for the 
original system. 


1/G,(1+8/G)) 


1 
c' = eafGses/GGssijG; *ite/G.40/G,G,t1/G,” (7) 
lt+s Gj ts G,G,t1 G, lt+s G,ts G,G, 1 G, 


Repeating equation 5: 
s/G,G, 


: 
© = 146/G,+8/6,G, 18/6, +4/G,G," 
l+s G,ts GG, l+s Gjts GG, (5) 


Comparison of these relations shows two 
differences: 


1 A term, 1/G,, has appeared in the 
denominator of both terms; if G5 is a high gain at 
tracking frequencies, this term will have little 


effect on the tracking response. 
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Bg The numerator of the disturbance 
term has been changed by the addition of the term 
1/G,; this gives a considerable reduction in the 
system isolation performance. 


Since the first effect is usually negligible 
or can be made so, the second effect is of 
principal interest. 


Figure 7 shows the response of both the 
original system (solid line) and the decoupled 
system (dashed line) to a severe transient to U. 
This figure shows the degradation in isolation 
performance which occurs. The transient used 
will be described under ''Practical Use of 
Decoupling. "' 


Review 


In this section, decoupling was explained 
by direct application to a typical system. The 
use of decoupling had two principal effects: 


Ife It eliminated the dynamic interaction 
between inner and outer loops, thus easing the 
stability problem considerably. 


Z. It reduced the isolation capability of 
the system by a considerable amount. 


In other words, the use of decoupling 
improved the stability characteristics at the 
expense of the isolation capability. In some 
applications, the loss of isolation capability is 
negligible. The next section describes a case 
(which uses the same system discussed here) in 
which the advantages of decoupling more than 
offset the reduced isolation performance. 


Practical Use of Decoupling 


System Configuration 


The system discussed in ''Application of 
Decoupling" (see figure 2) is typical of the type 
used for angle tracking in airborne radar. The 
previous section has shown the effects which 
decoupling will have on a system of this general 
type. This section indicates the effects which 
will occur when decoupling is applied to an 
airborne radar angle tracking loop of this same 
general type. 


In an airborne radar, the angular rate 
accuracy is usually more critical than the 
angular position accuracy. This fact, combined 
with the stability considerations covered 
previously, makes decoupling a useful technique 
in the design of this type of loop. 


The system shown previously as figure 2 is 
used in the following discussion. Its block 
diagram is redrawn as figure 8 to show the rate 
output. 


Original Loop Stability and Performance 


Since this loop is identical to that used in 
the discussion under ''Application of Decoupling, '"' 
the stability and performance characteristics are 
identical. However, the rate signal is now of 
primary interest. Equations 8 and 9 describe the 
position and rate performance. 


l SUS 
lt+s G,ts/G,G, l+s G,t8/G,G, 
s(1+1/G.,) s/G 
LS EEE aN = een ee Ty (9) 
1+8/G,+s/G)G, 1+8/G,+8/G,G, 


Equation 8 is identical to equation 5 which 
has been discussed in the section entitled 
"Application of Decoupling. "' 


Equation 9 which describes the rate 
performance of the loop shows several interesting 
features. 


Le The terms s/G> will tend to amplify 
high frequency noise in the input. 


(50 The disturbance present in the 
positional term appears in the rate term 
multiplied by G,. 


This indicates that the rate term is much 
more sensitive to both disturbance, U, and noise 
in the input signal than is the positional signal. 
Since this term is used (after multiplication by a 
prediction factor which is usually greater than 
unity) to predict target future position in airborne 
fire control systems, the presence of these noise 
signals is undesirable. 


These noise considerations, combined with 
the low stability of the system (using the same 
assumed transfer functions as before - see 
figure 4) indicate that this loop design is not 
satisfactory. In this instance, the use ofa 
decoupling loop will give a considerable improve- 
ment in system performance. 


De coupled Loop Stability 


Since the same configuration and transfer 
functions have been used here and in the section 
entitled ''Application of Decoupling,'' the same 
change in stability will occur. The addition of the 
decoupling loop will make the overall system act, 
dynamically, as if two independent loops were 
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operating in cascade. This is illustrated in 


figure 9. 


The only stability requirement for the 
decoupled system is that Gp and G)/s, separately, 
be stable open loop transfer functions. The two 
transfer functions no longer interact and may be 
designed separately. 


The transfer functions assumed here are 
such that the addition of the decoupling loop 
completely solves the stability problem. In 
addition, it permits a redesign of G; and G2 to 
give practically any desired characteristic. 
However, this paper will not explore the 
possibilities in redesign of the transfer functions. 


Decoupled Loop Performance 


The rate performance of the angle tracking 
system has been improved by the decoupling loop 
while the positional performance has been 
degraded. This is indicated by the response 
equations. 


1 


C!s R + 
(141/G,) (1+8/G,) 


(10) 


s 


vi "T¥s/G, * 


The effect of platform disturbance has 
been completely eliminated from the rate term 
but its effect on the position term has been 
increased by a large factor (see figures 10 and11). 


(11) 


Figures 7, 10, and 11 were obtained by 
analog simulation of a complex angle tracking 
control system. The disturbance was adjusted to 
be equivalent to a large rate disturbance. To be 
more precise, the disturbance used was 
equivalent to that seen by an airborne tracking 
radar with a lead angle of s 20 degrees, when the 
aircraft rolls at 100 degrees per second. 


Review 


"Practical Use of Decoupling" has 
described the application of the type of feedback 
compensation we have termed decoupling to a 
simple loop. The effects obtained were, of 
course, the same as were indicated in ''Appli- 
cation of Decoupling.'' However, several 
additional considerations became apparent. 


is The controlled variable was not of 
principal interest, and, although the position 
control was degraded, a complete cancellation of 
the disturbance effect on V was achieved. 


(ae Since the decoupling makes the loops 
act, dynamically, as if they were independent 
and operating in cascade, the two transfer 
functions can be designed independently. 


Bho The degradation in output positional 
accuracy is considerable and care must be taken 
to ensure that this does not become excessive. 


Component Matching 


General 


The preceding discussion is quite 
restrictive in that the transfer functions used 
have been oversimplified in order to describe the 
technique. Practically speaking, there must be 
transfer functions in both feedback paths of 
figure 5 to represent the physical elements 
necessary to convert mechanical motions to 
electrical signals, transmit the signals, and so 
forth. These elements are usually not 
represented in simple systems, since their 
effects can, generally, be absorbed in other 
transfer functions without loss of realism. 
However, in the case at hand, absorbing the 
effects of these elements in other transfer 
functions or omitting them from consideration 
implies severe restrictions on the physical 
elements involved in the feedback paths. That is, 
the extent to which a given loop can be decoupled 
and the utility of the resultant loop are limited by 
the characteristics of individual components and 


by the extent to which components can be matched. 


Feedback Transfer Functions 


Figure 5 can be redrawn, as illustrated in 
figure 12, to include the effects of the feedback 
transfer functions. We have shown the decoupling 
feedback as going to the internal summing point 
in order to represent more closely the mechani- 
zation which would actually be used in an airborne 
tracking loop. 


This tracking configuration gives the 
following response equations: 


G,G,G 1+G,G,G 


Le 395 
C= R+ we tl 2) 
+G G 
1G, G{G,G,46,G,G,G,G,G," 1+G,G,G,G,+G,G,G,+6,6,G, 
Ma G,(1+G,G,G,) ae G(G,G,G,) - GG, 4 (13) 
x + + +G,G_G 
1+G)GG,G,+G,G,G4tG,G,G,° 1+G)G)G,G,+G)G,G,+G,G,G, 
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These equations can be rewritten in a form 
suitable for direct comparison to equations 10 and 
11 as follows: 


1/G, [S,+ 1/G,G,] 


1 
c'= Gynl/GiGicmar/Gmar/ai + U (14 
Ge41/G,G,6 +6 ,/6,46,/G, © G,41/G,6,6,+G,/G,+G,/G, ) 
G e 
ee) 4t/G,G, % 1/G, [.S, G,/G,| 3 (15) 
G{H1/G,G,6+G,/6,4G,7G," * G76, G56,+6,/6 +G,/G, 


If equations 14 and 15 are compared to 
equations 10 and 11, several differences become 
apparent, and the condition which the individual 
transfer functions must meet, if successful 
decoupling is to be accomplished, is readily 
derived. 


Requirements for Decoupling 


Consider, first, the terms which provide 
the cancellation of disturbance, U, in the rate 
term, V. These terms appear as the numerator 
of the transfer function for U in equation 15. 


1/G, [<,S, : G,/G,| 


2 " GA/G,G,G,*G,/G+G Bo a 


( 

¥. 5 oS 

If V is to be unaffected by U, we have the 
requirement 


GiGae eG H/Ge (17) 


If the loop is to be decoupled, the 
characteristic function must be separable into two 
factors, one factor independent of G, and the 


other independent of G,. This will occur 
provided: 


G,G,G, = Ge 


(18) 

This is exactly the same requirement as 
was obtained for the elimination of U from the 
rate signal. 


The function, G,, is not available to the 
designer of an airborne angle tracking loop, since 
it is operating on the space positional feedback. 
For purposes of discussion, let us assume that 
G. is adjusted to equal the obtained value of Ge 
and that G, is adjusted to equal 1/G,. 

(Note that here we are artificially introducing a 
double matching problem; only one matching, 
given by equation 18, is basically required. ) 


1/G 
jae 1 2 
© "(GG G,) 141765) * "T1176 ,)" (19) 
1/G, 


G,+1/G,G, 


Vv! R (20) 


If these two requirements are imposed, the 
decoupled loop has the dynamic equivalent shown 
as figure 13. 


The equations from which the diagram is 
derived follow directly from the imposition of the 
above conditions on equations 14 and'15. 


Figure 14 illustrates the dynamic 


equivalent which is obtained if G,G,G is taken 


5 
equal to G; without the further assumption used 


to obtain figure 13. 


The equations from which figure 14 is 
derived are as follows: 


c (21) 


aa L R + cee U 
feat (Ce/ Sat Sy G/G,4/6, 


1/G 1/G,:O-U 
So dal Peale eae (22) 
G,+1/G,G, G,/G,+1/G, 


Review 


As indicated by equations 18 through 21, 
the decoupling of the loop is strongly affected by 
the extent to which certain transfer functions are 
matched. This matching problem is not, 
however, as severe for the decoupling approach 
as it normally is for matching filters, feed- 
forward, and conditional feedback techniques, 
since, in this case, the matching is done within 
the control loop, and the effect of mismatch is 
attenuated by the loop gain (see equations 14 and 
15). 


The loops as shown in figure 14 are not 
completely decoupled, since both loops involve 
the function Ge. Note, however, that the loops 
can still be designed independently provided that 
G, is matched to Gg. Even if these two 
components are not matched, considerably more 
latitude is available in the design of G, and G, 
than is available without the use of a decoupling 
loop*. 


*Throughout this paper, it is assumed that the 
transfer functions completely describe the loop 
characteristics, including loading effects. | 
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Component Drift 


General 


Precise noise and drift free components 
must be used in the decoupling feedback path. 
The following paragraphs discuss the effect of 
drift in various components and serve to provide 
a mathematical basis for this statement. 


Figure 15 illustrates the drifts to be 
discussed. 


Transfer Functions 


Direct solution of the block diagram shown 
as figure 15 gives the following response 
equations: 


“ 4 + [14+G,G 
G,G,G,R -G,G,G,G,D, -G,G,G,D, - G,G,G,G,D, fi ; 32s] U (23) 
ee + ‘ 
[1#G,G,G,G, +G,G,G, G,G,c,] (o} 
G,(1#G,G,G,) R - G,G,(14G,G,G,)D, + G,G,G,(G,+G,G,)D, 
é z 24 
G,G,(1+G,G,G,)D, + G\(G,G,G,-G,) U ( ) 
— 1 
=(14+G,G,G,G,+G,G,G, +G,G,G,) V 


If we now assume that the decoupling circuit 
is designed properly (i. e. , that G3zG4Gz = Gy), 
equations 23 and 24 reduce to the following form. 
G 


D (/G,) U 


1 51 Cpe Sele 25 
Crete he ee + (25) 
D D D D Gy. G,+1 G, 

Q/G,)R — (G,/G,)D, G,G,D, (G,/G,) D, (1/G,) ou 
ve ae es (26) 


h 2 7D-= 
where (Ge+1/G)G,)(G./G,+1/G,) 
DUES ely AG: 
5 1/ ts 


For comparison, the equations describing 
drift performance for the original loop are: 


C=G GjG_R+tU-G,G.G .D,-G.G\G,G,D 


[1#6,6,6,6,+6,6,c,] 
1 hae x 233) dee) tees Ons 


23) Gaoaes. + 


(27) 


G.G 


[I +G,G,G,c +G.G.G 
2534 


lee omcae3 4] ¥ = G, +e 


U+G,.G,G.G,G,D 


Rt= 
)R - G\G,U+G,G,G,G,G,D, 


28) 


- GG, [1 +6,c,c,] D, 


These can be put in a form more suitable 
for comparison to equations 25 and 26 as follows: 


63 (29) 
F 


where: F = (1/G,G,G,) + G,/G, + Gy 


Compa rison 


Table 1 summarizes the performance of 
both the decoupled and the original loops. 


Table 1 is not of immediate assistance in 
obtaining a comparison between the original and 
the decoupled loops, since, in the general case, 
a wide range of parameter values could be used 
for each of the transfer functions. If the follow- 
ing assumptions are made, a more obvious com- 
parison, which is typical of airborne angle 
tracking loops, is obtained. 


The result of applying these assumptions to 
the relations of table 1 is shown as table 3. 


Review 


For the angle tracking loop defined by the 
assumptions of table 2, the addition of the 
decoupling network leaves the drift character- 
istics essentially unchanged (G2 is generally a 
high gain), except that it provides an additional 
loop within which drift can occur. 


The effect of drift on the more general 
configuration depends upon the characteristics of 
the transfer functions in the manner shown in 
table 1. 


TABLE 1 


COMPARISON BETWEEN ORIGINAL AND DECOUPLED LOOP RESPONSE EQUATIONS - GENERAL 


C RESPONSE 
ORIGINAL 


* Fis 1/G,G,G, +G, + G,/G, 


6 


* 
% 
[ss 
il 


[¢. . 1/G,G,| [S6/Ss + 1/G,] 
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C'RESPONSE 
DECOUPLED 


V RESPONSE 
ORIGINAL 


V' RESPONSE 
DECOUPLED 


G.(1+G,G,G, ) 


Cee les 
sare 


is Wes 1 ste okeok 
3 D 
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Fig. 3. Decoupled angle track loop dynamic equivalent. 
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Fig. 4. Approximate gain-frequency curve for 


assumed tracking loop. 


Fig. 5. Angle tracking loop, decoupled. 


Fig. 9. Decoupled airborne angle tracking loop 
dynamic equivalent. 


Fig. 6. Decoupled system—dynamic equivalent. 
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Fig. 7. Response to transient disturbance. 


Fig. 10. Response of original system to a transient 
disturbance. 


Fig. 8. Airborne radar angle tracking loop. 
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Fig. 12. Angle tracking loop, decoupled. 


Fig. 13. Decoupled angle tracking dynamic equivalent. 


Fig. 11. Response of decoupled system to a transient 
disturbance. 


Fig. 15. Block diagram including drifts. 
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OPTIMUM COMPENSATION OF A POSITION SERVO 
WITH A MAGNETIC CLUICH ACTUATOR 


Ronald J. Hruby 
Northrop Corporation, Norair Division 
Hawthorne, California 


Summary 


In the past decade improvements have been 
made in the design of dry powder magnetic 
clutches which now make them competitive with 
hydraulic devices for missile servo actuator 
applications. Some of the inherent advantages of 
magnetic clutch actuators are: (a) increased 
control system reliability, (b) simplified equip- 
ment requirements, (c) properties remain constant 
with usage, (d) environmental requirements are 
quite compatible with an orbiting vehicle. The 
dynamic servo response of the magnetic clutch 
servo actuator is required in order to incorpo- 
rate these advantages in missile and space con- 
trol system design. 


This paper presents a synthesis procedure 
for the optimum design of a position servo based 
upon the closed-loop transient response to a step 
input. The transient response of the position 
servo is described by four independent para- 
meters. These parameters are the rise time, 
damping ratio, undamped natural frequency, and 
the steady-state gain. The synthesis procedure 
is based upon a theoretical model of a dry powder 
magnetic clutch, which is an air-gap, iron-care 
transformer with a one-turn secondary. The air 
gap is filled with dry ferrite particles. To 
provide simple analytic functions, a definition 
for rise time is developed which corresponds to 
the first crossing of unity gain to a high order 
of accuracy. The theoretical model of the mag- 
netic clutch has two, first-order time lags. 
These two time lags are derived from the proper- 
ties of the theoretical clutch model. They arise 
from the inductive properties of the clutch exci- 
tation coil and the clutch rotor induced eddy 
current, which is the current in the equivalent 
transformer single-turn secondary. The compen- 
sation for the excitation coil time lag consists 
of excitation current feedback, which decreases 
the effective transfer function time constant. 
The induced rotor eddy current time lag is opti- 
mally compensated by the proper selection of 
feedback gains as outlined in the synthesis 
procedure. 


The position servo drives an undamped iner- 
tial load. This is characteristic of a control 
system for an outer space vehicle or a guided 
missile which utilizes swivel rocket nozzles for 
control. It is also a first-order approximation 
to an airframe control servo utilizing neutral 
aerodynamic surfaces with low damping. The posi- 
tion servo characteristics are provided by a 
position feedback branch, and the basic servo 
stability is provided by a tachometer feedback. 
The theoretical model results in a third-order 
system which cannot be unstable. Practical mag- 
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netic clutches differ from the theoretical model 
primarily because they have ferromagnetic 
hysteresis. At high frequencies, hysteresis 
appears as a constant time delay which makes 
servo instability possible. An additional con- 
straint to the synthesis procedure is given which 
guarantees a stable, closed-loop position servo 
with the desired performance characteristics. 


The entire analysis in this paper is based 
upon the application of Root-Locus analytical 
methods to extract the theoretical roots of the 
third-order system describing the magnetic clutch 
actuator position servo. 


Introduction 


This paper presents a design philosophy 
which differs from the modern trend in servo 
analysis. A mathematical model of the position 
servo is formulated, which allows the exact 
synthesis of the closed-loop servo, based entire- 
ly upon the desired performance characteristics 
rather than the stability requirements of gain 
and phase margin. The mathematical model is 
based upon the simplification of the transfer 
function of the magnetic clutch by using an 
equivalent model of a transformer with a one-turn 
secondary. The result of this design procedure 
is an analytic solution of the third-order system 
for each of the system parameters, such as feed- 
back and forward branch gains as a function of 
the ideal performance characteristics. The final 
step in the design is a stability analysis. 


2 


Drivi tic Clutch Unit? 


lifiers and Ma 

Magnetic clutches used in missile servo 
actuators consist of an excitation coil which 
creates a changing magnetic field in a ferrite 
particle clutch gap. The existence of the mag- 
netic field allows the transfer of torque through 
the clutch to the output shaft. 


The basic electrical-mechanical configura- 
tion of the amplifier-magnetic clutch unit is 
shown in Figure 1. An a-c or d-c motor is the 
torque driving unit. It is mechanically coupled 
to the torque input spur gears on each clutch. 
The output torque shaft of each clutch is mechan- 
ically coupled to the servo output shaft. The 
clutches rotate in opposite directions to each 
other. Clutch #1 provides counterclockwise 
torque and clutch #2 provides clockwise torque. 
The two driving amplifiers each consist of a two- 
input difference amplifier driving a pentode 
class-B output stage. The output voltage for 
each amplifier varies from a small quiescent 
value to a maximum positive value. The input 
difference amplifiers are connected so that phase 


inversion is accomplished. When a positive volt- 
age is applied to the input of the unit, the cur- 
rent in clutch #2 approaches zero and the current 
in clutch #1 increases so that the output shaft 
has a counterclockwise torque. If a negative 
voltage is applied to the input of the unit, the 
current in clutch #1 approaches zero and an in- 
creasing current flows through clutch #2 which 
provides a clockwise torque to the output shaft. 
The amplifier clutch unit operates as a class-B 
push-pull system. With zero input voltage, a 
quiescent current is established in each clutch 
coil which semi-locks the output shaft position 
at zero output torque. The torque available at 
the output shaft is due to the drive motor. 
Changes in either the motor or load characteris- 
tics are in no way reflected back to the clutch 
currents via the clutch body. 


The voltage developed across the resistor 
(R) is fed back to the input difference amplifier 
for compensating the excitation coil time lag. 


The problem of inductive loading is severe 
enough that the use of pentodes to make a voltage- 
to-current amplifier without feedback in place of 
a voltage-to-voltage amplifier with feedback is 
not adequate because the self-inductance of the 
clutch coils will reduce the plate voltage until 
the pentodes operate in a quasi-triode mode dur- 
ing current transients. 


Analysis and Servo Co! nsation of Clutch Coil 


Figure 2a shows the equivalent clutch elec- 
trical circuit and equation 1 shows the resulting 
transfer function for (i/e}). 


Aig. oe ‘ 1) 
6 = RS rae : 


wv 


internal coil resistance 
K, = gain of voltage amplifier 
[. = excitation coil self-inductance. 


Tes te 


Ko = 


Figure 2b shows the method of decreasing the 
excitation coil time constant (“¥) to any degree 
by feeding back a voltage porportional to excita- 
tion current. The resulting transfer function 
relating the excitation current (4) to the input 
voltage ( @,) is given in equation (2). Both 
Figures 2a and 2b are equivalent circuits for the 
hardware system shown in Figure 1, without and 
with coil time lag compensation. 
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Rae natal apa 4 (2) 
~@ — Re RE+Kr 1+ RRIF Hl 


Spt 
where: 


ae as 
To * RF ROPER 
K, > RoFRE oR 


gain of current feedback network 


ae 


R = excitation current measuring resistor 


If the feedback term (Kgh) is large, the 
time constant (%) for the coil and amplifier is 
small. There is a restriction that R be small 
to minimize the input power requirements. By ad- 
justing Ko, h and R the value of (TQ) can be 
made as small as desired. The feedback network 
gain (fh) will be less than 1, and the network is 
required to provide a d-c path from the output 
side of (R) to the grid of the input difference 
amplifier, as well as to provide proper biasing 
of the input amplifier grid. 


Derivation of Rotor Eddy Current Time Lag 


The second first-order time delay (Te ) 
arises from the induced secondary current in the 
clutch rotor, which makes the entire unit act 
like a poor transformer. The rotor has conduct- 
ing elements by the nature of its structure, and 
the elements are equivalent to a single-turn 
secondary. The interaction of the driving cur- 
rent and induced secondary current produces the 
experimentally measured flux time lag. 


2 


Neglecting hysteresis, viscosity, bearing 
frictions, windage and drive motor transients, 
and assuming the drive motor is ideally com- 
pounded, the relationship between the magnetic 
flux and the exciting coil current is given by 
equations (3) and (4). It is a very good approx- 
imation to assume that the shaft output torque 
is directly proportional to the magnetic flux. 
This approximation is derived by neglecting the 
drive motor loading transients and using a push- 
pull clutch arrangement. 


Equations (3) and (4) are derived by treat- 
ing the rotor as a single-turn secondary with an 
equivalent winding resistance (Rs) and summing 
all the voltages and magnetomotive forces 
separately. 


Sum of all voltage sources in the equiv- 
alent secondary: 


. Ms (3) 
1sRs =e 10992 =) 


Sum of all magnetomotive forces in the 
netic circuit: 


mag- 


Nit ts Satis (4) 


where: (cgs - electromagnetic units - 
unrationalized) 
e 
1 - coil exciting current (amperes) 
ts = single-turn secondary current (amperes) 


single-turn secondary equivalent 
resistance (ohms) 


mean flux path length thru total clutch 
gap (centimeters) 


magnetic circuit mean cross-sectional 
area of clutch gap (sq. cm.) 


permeability of clutch gap 


number of primary turns (turns) 


flux to torque coefficient (ft. lbs./ 
Maxwell) 


Vv 


vy 


All magnetomotive forces are developed in 
the clutch gap, since the permeability of the 
ferrite core is much greater than the permeabil- 
ity of the clutch gap. The permeability of the 
rotor is taken to be the same as the gap energiz- 
ing ferrite core.? 


magnetic flux (Maxwells) 


The shaft output torque is given by: 
eat (5) 


Combining equations (3), (4) and (5) gives 
the results exhibited in equation (6), which shows 
the transfer function relating output torque to 
coil current (T/{), and the first-order flux time 


lag (Te )- 
By = 4 a2 Seer Ve 6 
Fea lego Oba A . 
& Rs 
ee 
— T+ “%P 


where: 


-8 

0). B fe) 

Te = LAS 

Ka SOS ee 

VAcuum: JL = 4,0 

Clutch Gap: JL~ SO To S00 
FERRITE CORE: U~10° To 10% 


PPM. 


The parameter (Vv) depends upon the driving 
motor torque characteristics and the clutch-motor 
gear linkage as well as the shear properties of 
the iron particles. (WV ) would have to be eval- 
uated experimentally. The flux time lag is pro- 
portional to the clutch area but inversely propor- 
tional to the rotor resistance ( Rs). 


The overall transfer function relating input 
voltage (©@,) to output torque T is given in 
equation (7). This transfer function includes 
compensation for ( Y ). 


aw Ka 


oF Cinctenase nds AA TILT ere? (7) 
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Optimum Synthesis Procedure 


In the preceding analysis the derivation and 
compensation of ( ‘YL ) was developed. Because 
the value of ( Tq) can be made very small by 
adjusting Ko, h and R, it will be set equal to 
zero for the closed-loop synthesis. This proce- 
dure is justified because the complex closed-loop 
roots are relatively unaffected by open-loop poles 
which are an order of magnitude greater than 


“). 


An optimum synthesis procedure must start 
with a set of performance parameters and provide 
analytic functions which compute the values of 
the various elements that constitute the closed- 
loop servo. ‘These performance parameters must be 
related to some type of exciting function which 
has been standardized by usage. The synthesis 
developed in this paper is based upon a unit step 
exciting function and four performance parameters 
which are the natural consequence of the unit 
step driving a finite servo. 


The damping ratio (©), the undamped natural 
(Wy), the rise time (ts), and the steady-state 
response (Eog,¢,) are the four performance para- 
meters used in the synthesis procedure. The 
damping ratio (5), the undamped natural frequency 

» and the steady-state response (Eog¢) are 
universally used performance parameters, but the 
rise time (ts) does not have a universally 
accepted definition. In this paper the rise time 
will be defined as the time required for the sin- 
isoidal term in the time domain to reach its 
first zero. For minimal values of ( ts), the 
output of the servo lies between 95% and 100% of 
its steady-state value at its first crossing for 
values of (§) between (.05) and (.95). Each of 
these performance parameters is shown perspec- 
tively in Figure 3. The definition for rise time 
(ts) was chosen because the synthesis equations 
are closed and analytic with respect to the 
design parameters. 


The complete closed-loop servo is composed 
of a forward branch which consists of an ampli- 
fier and clutch unit as shown in Figure 1, an 


inertial load (J), and both position and tach- 
ometer feedback branches which have gains ( ()) 
and (Q)), respectively. The complete servo is 
shown in Figure 4. The transfer function relat- 
ing the output position (€@,) to the input signal 
(@4) is given in equation (8), which includes 
the preceeding assumptions. 


eu Ki Ka (8) 
© TI Pp? VF uP KE Ki K203P + KikK2@y 
where: 
= load inertia (slug - ft2) 
Te = clutch time constant (1/sec) 
= amplifier-clutch coil gain (amps ) 
(volts) 
K> = magnetic clutch gain (ft 1bs/amp) 
Q,= position feedback gain = (volts/radian) 
Q,= tachometer feedback gain (volts/sec.x 
a 


radian) 


This transfer function can be resolved into 
three roots--two complex ones and one real root. 
The two complex roots are described by the tran- 
sient specification (¥ ) and (Wy). The real 
root is defined as (¥). The equivalent transfer 
function using (%), (3), (Mp) and the required 
steady-state value of (@,) (which was defined as 

©€occ), is given in equation (9). 


(9) 
Qe) - sex 8 wh 
Ci) ~ PP + +2 Foy’ + [2 S¥untrus]P+sw% 


The synthesis procedure can be performed in 
two ways. First, the four performance parameters 
can be specified and the values of (Qj), (Q)), 
( KsK2), and (%) can be computed, or a speci- 
fic clutch with fixed values of (‘Y¢) and (K,) 
can be used to compute therise time (ts), ( @,), 
(Az), and (K,) for specified values of (5) and 
(wy). If a specific clutch is used the compen- 
sation of (‘T._) is completed as the last step us- 
ing equation (2) and solving for ( Ke), using the 
criteria that[10T < WW]. ‘The exact value 
of (Yq) is immaterial provided [ma<tAo] - 
The synthesis procedure is derived by applying 
the Root-Locus gain and phase criteria to equa- 
tions (8) and (9) with the restriction that the 
complex roots of equation (8) be specified by the 
performance parameters (§) and (Wy). 


The relationship between the rise time (ts) 
and the real root (%) is one of the keys to the 
synthesis of the servo. This relationship is 
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given in equation (10). 
(10) 
OmniVi- $2! 
Tan(oVi- $?-t, +Tat fa ar } 


The real root (&) is restricted to be great- 
er than the real part of the complex roots which 
then provides a minimum and maximum obtainable 
value of ({g). The minimum value of (ts) occurs 
when ( ¥—> co ) and corresponds to a simple 
second-order system. The maximum value of (t5) 
corresponds to ( ¥—> Sm). ‘The original per- 
formance specifications will provide the maximum 
acceptable (ts) which in turn will specify the 
value of the real root (&%). There is a fixed 
relationship between the value of the closed-loop 
real root (6) and the magnetic clutch time con- 
stant (Yc). This relationship is given in equa- 
tion.(11):; 


8 = FSOn+ 


w= 4 (11) 
CT ¥FRSUm 


Equations (10) and (11) provide a means for com- 
puting the required (Te) given a specified rise 
time (ts). ‘The restrictions on (¥) place an 
upper and lower limit on (tg). The minimum and 
maximum values of (tgs) are given in equation (12) 


in terms of (§) and (Gy). 
(12) 
Usmin = 


Usimax = 


An alternate synthesis procedure is to start with 
a given clutch in which (Tj) is given. The 
value of (¥) can then be computed from equation 
(11), and the closed-loop servo rise time can be 
computed from equation (13). The real root (¥) 
is again constrained to be greater than ( SQ). 


te (13) 
te = Tea ES) pg 
S= GnvVi- s© 


With either method of fixing the value of (Te ) 
the synthesis can be completed by computing Q, , 
Gz and ( KifKa). An examination of equations 
10) and (11) immediately fixes the value of 
ey This is shown in equation (14). 


4 
Cogs, 


(14) 


Gas 


The application of the Root-Locus criteria 
leads to the evaluation of the tachometer feed- 
back gain (Qj). ‘The value for (Q@,) is given in 


equation (15) in terms of (7%), We (GEos¢ ), 
and (a). ‘The minimum value of (9, \ used in 
equation (15), 

(15) 


2 Soe ERED AN CRS See 
aoe Cogs (wrMi-3* + SenTan(d,) 


where: 


(2, = — Bes ts) Twi (HOVE) 


is determined by setting (t7=O). This corres- 
ponds to the second-order system with a minimum 
(te). With given or computed values of (%), 
($5, (FS), (Cy), (Gy), and (Wop), the value of 
( Ke) can be evaluated from the Root-Locus 
gain criteria. The value of the product ( KiK2) 
is given in equation (16). 


(16) 


[KiKs)] = on (1 —% SeonFeuke-$4))F 
2) ((Q,—02 Sum} + awh 1-82) 


If equation (16) is used to provide clutch 
specifications, then (7%) will be chosen and the 
value of (K,) can be computed by making suitable 
engineering choices for ({_) and (KK). If the 
choice for the clutch has been made, then ( Kz) 
is known and the required value of (Ki) is com- 
puted from equation (16). (Ky) is the closed- 
loop gain of the driving amplifier and the value 
of Ko) must be derived from the restrictions 
on (a.), R and fy for the fixed value of (K) ) 
and (Te ). Equation (17) shows this relation- 
ship. 


Kut (17) 


Ta 
7a. Ko => To SPECIFY CLUTCH 


Ko 
Ky 


=> CLUTCH GiveEN 


where: 
Ri, Tes L, Ka, AnD CK: Kel Known s 
CONSTRAINT 214 < 1/10 
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If specifications for a clutch are to be 
made, then (Te ) is computed from (tg) and the 
value of (K,) is determined from equation (17) 
by assuming nominal values of (Ko) and (L). 


The design of the excitation coil amplifier 
with a nominal gain (Ko) involves a locus of 
values for (R) and (h). The acceptable values 
of (R) and (hf) are determined by the value of 


(Ki) and (Ya). 


Amplitude Characteristics 
With Respect to Rise Time (ts ) 

The normalized amplitude of the transient 
response at the rise time (ts) is a function of 
the closed-loop parameters (5), (&) and (ts). 
For near optimum values of (3 ) the amplitude is 
near the first crossing of its steady-state 
value. Values of (8) less than (.5) would prob- 
ably not be chosen for normal servo synthesis, 
but the equations developed in the synthesis sec- 
tion are equally useful for analysis provided the 
exact relationship between the transient ampli- 
tude at the rise time (tg) for the prescribed 
range of (¥§) is known. The rise time (tg) has 
a minimum and maximum value for each pair of (3) 
and (GJy). By computing the transient amplitude 
as a function of a fixed percentage increase of 
(ts) over its minimum value, or a fixed percent- 
age decrease over its maximum value, the ampli- 
tude can be expressed independently of (Wm). 


Figure 5 shows the percent difference be- 
tween the amplitude of the transient and the 
steady-state value as a function of (ts/tsmin ); 
( ts/tsmax), and (3). Figure 6 shows the value 
of Pon as a function of (ts) and (¥). 
Figure 7 shows the relationship between 
( tsmin {OY ) and ( tsmex/Un ) for various 
values of - An examination of Figure 5 
shows that the percent difference in response 
amplitude is at least an order of magnitude 
smaller than the accuracy of the hardware com- 
ponents that would constitute the position servo. 


Servo Stabili $55 


Considerations 

In actual practice, the amplitude-phase 
characteristics of dry-powder magnetic clutches 
differ from a first-order lag because of magnetic 
hysteresis. 


Magnetic hysteresis appears as a fixed 
phase lag at very low frequencies and approxi- 
mates a constant time delay at high frequencies. 
The transfer function of the clutch at high fre- 
quencies has a phase lag which is porportional 
to the frequency. For these reasons the only 
complete approximation of the magnetic clutch is 
a meromorphic function. However, the amplitude 
characteristics at servo signal frequencies 
approximates a first-order lag very closely and 
the magnitude of the transfer function 
approaches ~6 to -8 db per octave at high 


frequencies. The main difference between the 
actual clutch characteristics and the equivalent 
first-order lag is the phase characteristic. The 
clutch model used in this synthesis procedure 
cannot produce an unstable servo. Since the 
physical clutch has a monotonic increasing phase 
lag, the position servo can be unstable. Figure 
8 shows a typical amplitude-phase plot of a dry- 
powder magnetic clutch. Indicated in Figure 8 is 
the first-order lag (Tc ) which would be used in 
the synthesis procedure. The first-order lag 
(Te ) is constructed from matching the amplitude 
of clutch transfer function to a theoretical lag 
without considering the phase plot. The equiv- 
alent first-order lag (7) should have a slight- 
ly higher or equal amplitude at signal frequen- 
cies than the clutch. The phase lag of the 
actual clutch is less than the equivalent first- 
order lag (Ye ) for all frequencies less than 
(Gq). At (GJ =Wa) the phase lag of the clutch 
and the equivalent lag (%%¢) are equal. For 
typical clutches (Wa) is greater than or equal 
to ( 3/%e ). These relationships are shown in 
Figure 8 


The synthesis procedure developed in this 
paper will always provide a stable closed-loop 
servo with actual clutches, provided the value of 
(Gm) is constrained to be less than or equal to 


(Wa), and (ts) is kept smaller than (Atginay)- 
These relationships are placed in their proper 
perspective when it is realized that (tg) would 
be taken as near ( tsgyin) 28 epagibiaie and 
(Gy) will usually be Yess than (1fzUX), the 
equivalent corner frequency of the clutch. 


References 


1. Magnetic Servo Clutch, February, 1951, Boeing 
Airplane Company, H. Lusissar and E. N. 
Zangen 


2. Characteristics of Some Magnetic-Fluid Clutch 
Servomechanisms, 1949, MIT Servomechanisms 
Laboratory, A. J. Parziale and P. D. Tilton 


3. Investigation of Some Design Fundamentals of 


the Magnetic Particle Clutch, July, 1950, 
MIT Servomechanism Laboratory, J. P. Ivaska, 


R. Kramer, and G. C. Newton, Jr. 


4, The Theory and Behavior of the Lear Magnetic 
-Powder Clutch Servo, September, 


1957, Lear Incorporated, A. Bepristis 


5. Basic Feedback Control System Design, 1958, 
McGraw-Hill, C. J. Savant, Jr. 


CLUTCH *! 


CLUTCH *2 


Fig. 1. Electrical-mechanical configuration of ampli- 


fier and clutch unit. 
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Fig. 2. (a) Excitation coil equivalent electrical circuit. 
(b) Equivalent circuit of Fig. 1 including coil 
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Fig. 3. Closed loop performance specifications. 
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Fig. 4. Complete closed loop servo with magnetic 
clutch actuator. 
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Fig. 5. Transient response at rise time (t,) with respect to the steady-state value Cog s.- 
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Fig. 6. Normalized real root (7) with respect to rise time (tg) and damping ratio (). 
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Fig. 7. Normalized minimum and maximum rise time vs damping ratio. 


S 
<< 
(O\2 
‘ So) 
ee oe 
pt oe 


TRANSFER FUNCTION OF TYPICAL 
MAGNETIC CLUTCH ACTUATOR 


30 -150 


TTP: pee te eis 
RS AMPLITUDE |_| 
PHASE 
RP 
TRANSFER FUNCTION OF EQUIVALENT 1ST PHASE 
ORDER LAG FOR MAGNETIC CLUTCH ACTUATOR | AMPLITUDE 
OSES ae ee 
NOTE: TO GUARANTEE SERVO STABILITY: 
ee aee a 
pL LL a 
c 
Ee 3 eelbas 5 ea 7s FO 2.0 30 40 5.06... 10, 


Fig. 8. Typical magnetic clutch actuator transfer function and equivalent theoretical first-order model. 
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Abstract 


The design of a conventional autopilot re- 
quires a detailed knowledge of the structural 
dynamics of the missile that is not easily obtain- 
able. The self-adaptive autopilot considered in 
this paper would not depend on this detailed know- 
ledge. This system would, ideally, adjust itself 
to maintain optimum performance in an environment 
of changing inertial, aerodynamic, and structural 
parameters. The mechanization of a suitable self 
adaptive system is complicated by the difficulty 
of getting a continuous measure of system perform- 
ance. The method used, in this study, to obtain 
the basic modes is based on decomposition of the 
transient response. 


Introduction 


In an autopilot system for a large flexible 
booster with relatively low bending frequencies 
and small structural damping the system parame- 
ters, such as rigid body natural frequency and 
damping, should be invariant to changes in the 
mass and aerodynamic environment of the vehicle. 
The design of a conventional autopilot involves 
a detailed knowledge of the structural dynamics 
of the missile that is not easy to obtain. There- 
fore an autopilot system must be evolved that 
does not depend on a detailed knowledge of the 
missile structural characteristics. 


The concept of the self adaptive autopilot 
has been applied in an attempt to solve these 
problems. Ideally, this system automatically 
adjusts itself to maintain optimum performance in 
an environment of changing inertial, aerodynamic, 
and structural parameters. The mechanization of 
a suitable self adaptive system is complicated by 
difficulties associated with the continuous meas- 
ure of system performance. A method based on de- 
composition of the transient response to obtain 
the basic modes is described. 


Satisfactory self adaption can be accomplish- 
ed if information on performance can be obtained 
with negligible time delay. This requirement ne- 
cessitated the development of a method for obtain- 
ing system performance parameters instantaneously. 
This method is based on the transient response of 
the system to periodic impulse disturbances. In 
general this response will consist of the rigid 
body modes plus the lightly damped elastic modes. 
The concept of the self tuning frequency sensor 
has been evolved to accomplish the decomposition. 
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This device is excited by the missile's impulse 
response and continuously furnishes information 
about the rigid body natural frequency and damp- 
ing. The frequencies of the elastic modes are 
also continuously furnished. Using this informa- 
tion, the rigid body dynamics of the missile are 
self adapted to be invariant with changing inertia 
and aerodynamic environment. The frequency in- 
formation on the elastic modes is used to adjust 
a compensating filter that keeps the structural 
feedback in the autopilot to a negligible level. 


This system can satisfactorily eliminate the 
structural feedback problem in a missile system 
including two elastic bending modes. The ability 
of the system to handle two bending modes having 
been demonstrated, the system can theoretically 
be extended to handle any number of bending modes. 


The Elastic Missile System 


Consider the block diagram of the general 
elastic missile shown in Figure 1, where 


plies +1 
. Diese ker 


n 
E tt : 
k=1 Ss + bs + 1 


a 1 
The closed loop transfer function for this system 
is 
AG 

E 


7) 9 
+ AG, (Kp ale Ks) 


Ss 


Substituting for Gas 


n 18° + Bs ae 


A I 7) 
8, k=1 as + bs +1 
85 : 2 n 1S pio at iL 
GrulgataDics tanta eae oer getapts 
k k 
Letting n = 2 and plotting the amplitude ratio 


and phase of this equation for typical values of 
the damping and natural frequencies associated 
with the elastic missile, it becomes evident that 
any theory for self-adaption which takes care of 
the elastic modes must be capable of distinguish- 


ing between these modes (Fig. 2). For example,if 
a self adapted, filter were successfully set up to 
cancel out the first bending mode characteristic, 
the problem would only have been transferred to 
the second bending mode, and so forth. Thus, any 
solution to the self adaption of the elastic mis- 
sile must be capable of taking care of several 
modes simultaneously. 


The self adaptive system must be able to: 


i) 
modes; 


Distinguish between the various elastic 


2) Adjust a filter for cancellation of each 
individual mode. 


To keep the problem within bounds, it was de- 
cided to study a rigid missile plus two bending 
modes. This is sufficient to demonstrate a method 
applicable for the general case where the system 
must: 


1) Distinguish between the first and second 
bending modes; 


2) Adjust a suitable filter for cancellation 
of these modes. 


Philosophy of Overall System 


The general approach to the solution of the 
problem is as follows: 


The rigid body mode is made self adaptive. 
The purpose of this feature is to establish the 
rigid body dynamics equal to some predetermined 
invariant factor. This provides the missile with 
uniform dynamics over the flight trajectory. The 
missile impulse response is then passed through a 
rejection filter that rejects the rigid body mode. 
This leaves a residue that consists of a super- 
position of the elastic mode responses. Since 
these modes possess essentially zero damping, the 
elastic residue can be expressed as 


n 


DL A. sin C(t ty). 


k=1 


fe) = 


The elastic residue is decomposed to obtain 
the component frequencies Ne The purpose of this 
operation is to obtain signals for the positioning 
of cancellation filters. Thus, since the elastic 
dynamics enter into the problem according to 


ys + B,8 + 1 
 _ 
2 k=1 as” + b8 +1 


a cancellation filter is reasonably employed that 
is essentially the reciprocal of this quantity. 
Examining the individual contributions from Gp 


(Fig. 2) shows that each mode has the effect of 
contributing a tall spike at the frequency Me 


Thus the approach to cancellation is based on 
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measuring the frequency and positioning a re- 
jeetion filter at the frequency ae This can be 
done having only a knowledge of Me 


In the determination of the elastic mode 
frequencies, the concept of the self tuning fre- 
quency sensor has been evolved. This device can 
be started with arbitrarily set values. It can 
tune itself by information received from the elas- 
tic residue to give out signals proportional to 
the component frequencies r in the elastic resi- 


due. If the problem requires knowledge of n bend- 
ing frequencies, n parallel channels of the device 
are employed. The theory presented is for two 
channels, but the device can easily be extended 
to more. 


The operation of the proposed autopilot sys- 
tem for two elastic bending modes is illustrated 
by the block diagram in Figure 3. As illustrated 
in this figure, the impulse response of the mis- 
Sile is applied to the self tuning frequency sen- 
sor. Since the rigid body mode is the lowest fre- 
quency, the first channel of this device selects 
and measures the rigid body natural frequency. A 
separate measurement of damping ratio is made, 
the theory of which is developed later.The values 
obtained for these quantities are then compared 
with the desired values and error signals 
obtained. These error signals are used to close 
a loop around the autopilot position and rate 
gains K, and K, to maintain the rigid body natur- 
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al frequency and damping ratio invariant. After 
rejection of the self adapted rigid body mode the 
remaining output signal represents the impulse 
response resulting from the first and second bend- 
ing modes, the elastic residue. This elastic re- 
sidue is the input to the second channel of the 
bending mode frequency sensor. The output from 
this channel is the square of the frequency of the 
first bending mode rejection filter, positioning 
this filter to suppress the first bending mode dy- 
namics. Simultaneously the first channel output 
positions a rejection filter for cancellation of 
the first bending mode from the input to the sec- 
ond channel of the frequency sensor. This channel 
works in exactly the same manner as the first 
channel, and selects the second bending mode fre- 
quency. 


The Frequency Sensor 


The differential equation describing harmon- 
ic motion is 


5g 2 
a+re= 0. 


The solution of this differential equation is a 
sine wave 


6n= Ae Sine tes 


‘ 


Thus, a sine wave and the second derivative of 

this wave are related by the square of the fre- 
quency. This principle can be used as a basis 

for measuring the frequency of an unknown sine 

wave. 


In the case of the elastic missile, the tran- 
sient response after removing the rigid body mo- 
tions is the sum of sinusoids of unknown frequen- 
cy and amplitude. A device is desired that is 
able to sort out these sinusoids and give informa- 
tion about the frequencies contained in the sum. 
The problem is solved under the additional restric- 
tion of allowing no prior knowledge of the compo- 
nent frequencies. 


Suppose the sum of sinusoids is run through 
a device whose output would be predominantly a 
Sinusoid of the lowest frequency contained in the 
sum, no matter what this frequency may be. The 
output of this device is then differentiated twice 
and the result of the double differentiation com- 
pared with the original. The magnitude of the two 
Signals would be in the ratio of the frequency 


squared. Suppose further that the signal nN is 


used to tune the original input filter to further 
favor the lowest frequency. In the limit, the fre- 
quency of the lowest frequency sinusoid in the sum 


2 
can be measured. Having a measure of Me this 


measure can be used to tune a rejection filter so, 
if the original signal is run through this filter, 
the lowest frequency component will be knocked out. 
If this signal is fed to a second channel similar 

to the first one described, it will select the si- 
nusoid of next highest frequency, Noo to favor and, 


in the limit, will measure Ne This procedure can 


be carried on as far as required. 


For details of the device, first consider the 
input filter that favors the lowest harmonic pre- 
sent. If this filter is made a simple undamped 
second order system, it has the transfer function 


k 
2 
ot i+ 4s 
2 ; 
where T is adjustable by Se The amplitude char- 


acteristic for this filter is the second order 
characteristic with zero damping (Fig. 4). 


Suppose initially 1/t is set to some value 
lower than any expected frequency present in the 
sum of sinusoids. Then the portion of the curve 
considered is on the right half of the infinite 
spike. This portion of the curve is asymptotic to 
-12 db/foctave. This will provide a two to one 
separation of amplitudes if the first and second 
frequencies are of equal magnitude and are separa- 
ted by two to one. Since, in the missile, the mag- 
nitude of the first bending will exceed that of the 
second, the output of the filter will favor the 
lower frequency. 


Assuming that the frequency can be measured, 


the filter can be tuned to further favor the lower 
frequency. If this process is carried to the limit 
for the system shown, infinite response to the low- 
est mode will result. Since this is undesirable, 
automatic gain control is used to limit the ampli- 
tude of this oscillation to some predetermined le- 
vel. This is done by feeding the output of the 
filter into a deadspace device, the width of which 
is equal to the double amplitude of the desired 
oscillation. If the oscillation exceeds the width 
of the deadspace, pulses will come out of the dead- 
space that can be fed back to limit the amplitude 
of the oscillation out of the filter. This intro- 
duces some distortion of the sine wave, but exper- 
ience has not shown this to be serious. 


The measurement of the frequency of the ampli- 
tude limited oscillation will now be studied. Con- 
sider the sine wave 96= A sin At. Differentiating 


this function twice 6 = ae sin \t. 


2 
Then =~ = -)\ . To accomplish the double differen- 


tiation, a network is used with a transfer function 
of 

Z 
s 


i) 
2 

i ie oie 

The lag term in the denominator is used to 
avoid having infinite gain at high frequencies, 
which would result in instability of the circuit. 
To provide the proper phase relations during the 
comparison process, the direct signal is passed 
through a network whose transfer function is 


al 
(le Se 11s)" 


If the direct signal were now to be multiplied 
by A and summed with the double differentiated 
Signal the result would be zero. Suppose, however, 


22 
that it is multiplied by some quantity, \' , that 


in general is not equal to \. The sum of the 
direct signal and the double differentiated signal 
is then a sine wave whose magnitude is equal to the 


errore = nr? - tt this is integrated with re- 

spect to time the result would again be zero.There- 
fore, the two signals are rectified before they are 
compared. This then gives the absolute value of 

the two waves, and on comparison a string of pulses 
is found. The pulses are positive or negative as € 
is positive or negative. If these pulses are inte- 


2 

grated and the output of the integrator oped ss 
2 

the stable condition for this integrator is \' =) , 


or ¢ = 0. This gives the quantity for adjustment 
of the input filter and to multiply the direct sig- 
nal by. A block diagram of one channel of the fre- 
quency sensor is shown in Figure 5. 


The Self Adapted Autopilot 


Consider now an autopilot system with provi- 
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sions for sensing and cancelling out two elastic 

bending modes. Such a system with no compensation 
for the elastic properties would be represented by 
the block diagram of Figure 2 and n = 2 in the term 


Gp Examining Figure 2, the open loop frequency 


response curves for one environment for this sys- 
tem corresponding to a typical large elastic boost- 
er, this system is seen to have a negative phase 
margin when the transfer function crosses the O db 
axis, and hence is unstable. The normal approach 
to this problem is to gather data allowing examina- 
tion of the stability of this system under all anti- 
cipated environmental conditions of the missile, 
and then design a compensating network to achieve 
stability. This procedure is objectionable from 
several points of view: 


1) It requires detailed knowledge of the 
anticipated missile structural and environmental 
characteristics. This is at best a very dubious 
and laborious task. 


2) It requires a compromise in the compensa- 
ting network that attempts to stabilize the loop 
over the expected structural and environmental 
ranges. This at best provides marginal perform- 
ance. The compensating network that will make the 
loop stable over the trajectory the missile is go- 
ing to fly is the result of a calculated guess. It 
is impossible to optimize the performance under 
these conditions. 


These objectionable problems can now be eli- 
minated. The first can be eliminated if the elas- 
tic modes contained in the missile response are 
sensed and this information used to adjust compen- 
sation to provide stable operation of the loop for 
any and all environments. If the elastic modes are 
effectively eliminated and nullified, only the ri- 
gid body response is left. If the rigid body por- 
tion of the autopilot is now made self adaptive, 
optimum performance of the missile response can be 
provided that is invariant to environmental condi- 
tions. 


As pointed out previously the elastic modes 
enter into the loop through a transfer function of 
the form 

2 
n 118 + Bs + 1 

CG. Se eee a edeinn ar , 
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The constants in this product are so related 
that each term contributes a tall spike in the open 
loop transfer function, the net contribution in 
phase being zero. Figure 6 shows a typical term. 


The altitude of the peak is of the order of 
+28db, due to the structural damping ratio of 0.02 
or less. If the location of the spike were known 
compensation could be used to provide cancellation, 
that is essentially the reciprocal of the mode 
transfer functions. Thus the approach is to sense 
the frequency of the bending modes and position a 
rejection filter in the loop to cancel the elastic 
structural feedback. Such a filter is shown in 


Figure 7. The success of the whole scheme thus 
depends on sensing the elastic frequencies Ne 


A block diagram of a self adapted autopilot 
for a large elastic booster, using the self tunin 
frequency sensor previously described, is shown i 
Figure 3. This figure shows a system for sensing 
and eliminating two elastic bending modes. 


So far the self adaption of the rigid body 
mode has not been described. This is done by man 
pulating the rigid body rate and displacement K, 
and Ky 
ing considered, the response characteristics of th 
inner loops (such as the hydraulic loop) and the 
dynamic characteristics of the sensing instruments 
is very high as compared to the inertial response 
characteristics of the missile. Thus, these re- 
sponse characteristics are separated by a factoro 
ten to one. Under these conditions, the dynamics 
of the inner loop and sensing instruments can be 
ignored and the assumption made that these element 
are essentially ideal. This leaves only inertia 
and aerodynamic forces as far as the rigid body is 
concerned. Under these assumptions the rigid tran 
fer function becomes 

2 


WW) 
G n 


respectively. In the class of vehicle be- 


; 2 ye 
i 5) AS25609 sa) 
n n 


If the elastic effects have been sensed and 
nullified, this becomes the overall closed loop 
transfer function of the missile autopilot system. 
In a space environment, aerodynamic effects become 
negligible and the damping and restoring forces on 
the missile are determined by the rate and dis- 
placement gains. To self adapt this transfer func- 
tion, the frequency wis sensed in a manner simi- 


lar to that described for the frequency sensor.This 
measured frequency is then compared with some de- 
sired set value for this frequency and an error 3 


n 
obtained. This error is fed into an integrator 


whose output is Kp: Thus, the integrator either 
increases or decreases Ky until w equals its de- 
n 


sired set value. A simple procedure can be follow- 
ed for the adjustment of the rigid body damping 
ratio. Consider the differential equation 


2 2 
oot 2tores er we Ts bes. 


Now suppose €, is the value of the damping ratio 


desired, and let& = &,+ At. Then 
Gottch (EF Ab)w 6) +)y 79 yo 
O ) BOs ee a ated oe 
Solving for w 0 AE, we have, 


2 nats 
w OAs ay (O.-G0) -2E Go. 


Since wis fixed at its known value, by mani- 


pulating Kp, and hence is known; theoretically At 
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could be calculated, if 6 ., @., P) ,» and 6 were 
(oy) ake Pos fe) 


known, by performing a division. This seems ob- 
jectionable, because it involves measuring four 
quantities and performing a division. What is de- 
sired is to measure two quantities and not perform 
a et This can be done as indicated in Fig- 
ure 0. 


The missile output or is measured and 8, and 


o. are generated by differentiation. 


is introduced in the differentiations to cut off 
the high frequency and provide stable operation. 
The output of the summing junction is a quantity, 
w 6 As 


A lag term 


that can vary in sign, not only with Aé 
2 
(l+ts) 


but also with 9,. Hence, this quantity is multi- 
plied by the output from the first differentiator 
to give a term whose sign varies only as At. The 
multiplier provides a quantity whose magnitude and 
sense is a measure of A€. If this quantity is 
integrated and the output of the integrator made 
Kp then this integrator either increases or de- 


creases Kp woegeat il 


At = 0. 


Thus, the autopilot system is self adapted to 
have dynamic response characteristics that are in- 
variant to the missile environmental conditions. 


Fig. 1. Elastic missile autopilot block diagram. 
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Fig. 3. Block diagram of self adapted autopilot for elastic missile. 
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Fig. 4. Input filter characteristic. 
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Fig. 6. Amplitude ratio for one elastic mode. Fig. 7. Rejection filter characteristic. 


Fig. 8. Damping ratio error detector. 
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DESIGN OF OPTIMUM BEAM FLEXURAL DAMPING IN A MISSILE 
BY APPLICATION OF ROOT-LOCUS TECHNIQUES 


Ronald J. Hruby 
Northrop Corporation, Norair Division 
Hawthorne, California 


Summary 


The control of a missile in free-fall in the 
presence of external disturbances requires that 
the control system be stable when various body 
bending modes are excited. It is an inevitable 
consequence of the internal and external missile 
environment that all body bending modes are 
excited. The general method of stabilizing a 
missile is to choose a minimum of inertial sen- 
sors such as rate and attitude gyros, and through 
the use of passive electrical networks in the in- 
formation paths, generate stabilizing poles and 
zeros in the control system open-loop transfer 
functions. 


Generally, the electrical shaping networks 
are approximations to "maximally" flat functions 
with the cutoff frequency below the first bending 
mode resonant frequency, and in extreme cases, a 
notch filter is needed to eliminate body bending 
information. The basic purpose of the electrical 
shaping networks is to attenuate the high signal 
frequencies circulating in the control loop with- 
out the loss of important low frequency data. In 
some applications the use of passive shaping net- 
works becomes very difficult because of the large 
attenuation required near the first body bending 
resonant frequency and the required amplitude 
flatness at lower frequencies. If passive shap- 
ing cannot be used, then active networks must be 
substituted. This is an undesirable complexity. 


This paper presents an optimum design pro- 
cedure to utilize the characteristic properties 
of completed airframe and missile structures to 
acquire a maximum increase in modal damping with- 
out auxiliary equipment and further autopilot 
complexity by the proper choice of the sub- 
assembly attaching fixtures. 


Completed airframes and missiles are com- 
posed of a very lightweight integral structure in 
which a large number of equipment sub-assemblies 
are attached. The sub-assemblies or equipment 
modules are generally massive and of small volume 
so that they act as a rigid mass attached to the 
lightweight integral structure via an equivalent 
lateral spring constant and viscous damper. For 
example, the guidance package and the rocket 
engine of a ballistic missile are usually attach- 
ed fore and aft of the missile integral tank 
structure. Each of these elements can be quite 
rigid and will act as a separate mass from the 
integral tank structure. The proper choice of 
the equivalent lateral spring constant and vis- 
cous damper of the attaching fixtures that hold 
the sub-assemblies to the basic lightweight mis- 
sile or airframe integral structure will result 
in an optimum increase in modal damping for a 
particular beam-flexural mode of the completed 
structure. ; 


237 


To determine the values of the proper equiv- 
alent spring constant and viscous coefficient, the 
dynamic principles underlying the Frahm damper are 
used by extending the analysis of the Frahm damper 
concept to include a general linear flexural beam 
which satisfies a linear partial differential 
equation which can be solved by the product of the 
solutions of two ordinary, linear differential 
equations. 


The extension of these principles is accomp- 
lished by deriving the solution of the attached 
sub-assembly to the integral structure through 
normal servomechanism analysis. This analysis 
uses the analytical relationships of the Evans or 
Root-Locus method. 


Introduction 


The design of air-space and outer-space 
vehicles usually results in light, rigid struc- 
tures which have the dynamic properties of beams 
with a high finesse ratio. The combination of 
rigidity, minimum weight, and slenderness means 
the first body bending mode, and possibly some of 
the other modes will have a very low damping ratio 
which will be on the order of 10% or less. 


The controls engineer is faced with the prob- 
lem of attaching inertial elements to the vehicle 
to provide flight attitude stability. Typical 
propulsion units such as rocket and jet engines 
that are used to drive air-space or outer-space 
vehicles also excite the structural flexural 
modes. 


The controls engineer must design his closed- 
loop control system so that the vibrational energy 
in the body-bending modes is not increased as the 
result of closing the attitude control loop. In 
effect, the control system must contain signal 
attenuation at body-bending frequencies so that 
the open-loop gain is less than one. The requir- 
ed amount of signal attenuation at body-bending 
frequencies increases very rapidly as the flex- 
ural mode damping ratio approaches zero. These 
requirements have led to overly complex autopilot 
systems with less than the desired reliability. 


This paper presents a passive method by 
which the damping ratio for any beam flexural 
mode can be optimally increased by taking advan- 
tage of the sub-system equipment modules which 
are attached to the missile body. Ina missile, 
the rocket engine assembly, the heavy payload 
package, or various electronic and guidance pack- 
ages can be used to improve the damping ratio of 
any of the body-flexural modes by choosing the 
attaching structural characteristics by the pro- 
cedure outlined in this paper. The analysis 
applies only to lateral bending modes since the 
longitudinal modes are relatively unimportant. 


The procedure for improving modal damping is based 
upon an extension of the Frahm damper principle 

to that of a free-fall beam with concentrated 
masses attached to the beam with an equivalent 
lateral spring and dashpot. 


Since the concentrated mass is necessary for 
the system operation and the spring is required 
for structural reasons, the dashpot is the only 
new element added to the missile system. The 
size of the equivalent dashpot is controlled by 
the concentrated mass and the modal resonant fre- 
quency. The weight venalty of the dashpot 
appears to be a negligible minimun. 


Description of a Frahm Damper System 


There are many mechanical structures which 
appear to be a large mass located at the end of a 
firmly imbedded spring. Such a structure is 
shown in Figure l. 


If a forcing function is applied to the sus- 
pended mass with a forcing frequency near the 
spring-mass resonant frequency, the resulting 
position amplitudes become extremely large and 
therefore undesirable. Theoretically, a forcing 
frequency equal to the spring-mass resonant fre- 
quency results in infinite position amplitudes. 
In a practical case, the mass is constrained by 
stops. The forcing function is usually an exter- 
nal disturbance on which no control can be exert- 
ed. The problem is to find a method of linearly 
restricting the positional amplitude of the mass 
without adding a damping device in parallel with 
the spring. The solution of this problem is the 
addition of a small mass, coupled with a spring 
and a dashpot to the large mass, with a natural 
frequency very close but not equal to the natural 
frequency of the large mass and spring. Such a 
system is shown in Figure 1 and is called a Frahm 
damper . 


The analysis of the simple system in Figure 
1 leads to equation (1), which gives the damping 
ratio (4m) of the system mass (M) in terms of 
the system parameters (K), (B), (K), and (m). 


Sa beaital* ~ 


zZ 
Mm +M. 


= Damping ratio of mass (M) 


ni- 


I 


= Spring constant of mass (M) (1bs/ft) 


= Frahm mass (slugs) 


238 


= Frahm spring constant (lbs/ft) 


K 
B 


The Frahm mass (™M) is arbitrary. In gener- 
al, the larger the Frahm mass (™M) the greater 
the increase in the effective damping of the large 
mass (M), but the Frahm mass (™) is usually 
restricted by other known considerations to be 
much smaller than (M). 


= Frahm dashpot coefficient (lbs. sec/ft) 


Application to a Free-Beam 


The usefulness of a Frahm damper to increase 
the modal damping of a free-beam and therefore a 
missile in free-fall, arises from the flexural 
characteristics of a beam which makes the differ- 
ential equations of motion similar to the system 
shown in Figure 1. Figure 2 shows how a Frahm 
damper is applied to a free-free beam. The damper 
is applied at a point (231) on the beam and is 
"tuned" to the mode of the beam which requires the 
increase in modal damping. The application of a 
Frahm damper to a free-free beam is like a servo 
feedback system in which the feedback gain is pro- 
portional to the normal mode amplitude squared and 
the mass of the Frahm damper. Since weight is 
precious for a missile, it requires the damper to 
be located at the point of maximum modal ampli- 
tude, usually at the tail or the nose (either end 
of the beam). ‘The analysis for a beam does not 
result in a set of simple equations such as exhib- 
ited in equation 1. However, a set of simultan- 
eous equations do result, which when solved give 
the same information required for designing the 
damper. The analysis of the problem is based upon 
a@ non-homogenous beam with low modal damping. ‘The 
solution is derived from knowledge of the beam's 
modal frequency (Wm), inherent damping (B,), 
orthogonal mode shape (GY(x)), modal mass (Mn), 
and the mass of the Frahm damper (m). ‘The anal- 
ysis is based upon the application of the Root- 
Locus theorems. 


Mathematical Analysis 


The basic differential equation of motion of 
a free-free beam is given by equation 2. This 
equation applies to the system shown in Figure 2. 
The solution of this equation depends upon the 
capability of forming a product solution of time 
and lateral displacement. 


The solution in time is taken to be a sinu- 
soid at the driving frequency of the forcing func- 
tion (i.e., the forcing function is taken to be 
the product of a sinusoid in time with a variable 
magnitude in (X)). This beam loading is illus- 
trated in Figure 3. 


Equation 2 applies to a non-homogeneous beam, 
bending in one plane with no gyroscopic or cross- 
dimensional coupling effects but including inher- 
ent beam damping. 


SE Oooo Rares) + (XB: 3609 ‘ 
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where: 

Foc) = Young's modulus (1bs/ft) 

Tex) = Cross-sectional moment of inertia 
(ft)4 = Res 

R =» Cross-sectional radius of gyration 
(£t)® 

S = Cross-sectional area (ft)2 

Bi = Damping force/mass/unit length 
(lbs.sec.ft./slug ft.) 

Aix) = mass/unit length (slugs/ft) 

fost = External force/mass/unit length 


(lbs. ft/slug) 


The lateral displacement y(x,t) is a product 
of a function of time and a function of (x). ‘The 
forcing function is also the product of two in- 
dependent functions. Equation (3) shows this 
relationship. 


MX,t) = YO 4 (4) (3) 


fex,t) = flac -$cp 


The function y1(x), in general, is a sum of 
independent functions of (x) each multiplied by 
an arbitrary constant which is evaluated by in- 
corporating the boundary conditions of the beam 
and the specific forcing function f}(x). Equa- 
tion (4) gives the general form yj(x). 


= (4) 
x) XL Gn Yn (x) 


where: 
- A modal constant which is derived from 


gn ficx) » YO and JMy(modal mass) 
Yp)= Orthogonal N-mode shape 
(00 is evaluated either by the analytical 


solution of the beam differential equation when 
possible or by experimentally recording the mode 
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shape which is the most generally used method. 

The analytical expression for (U(x) has as many 
unknowns as the independent boundary conditions of 
a free-free beam. These boundary conditions are 
that the third and fourth derivatives of the lat- 
eral displacement at both ends of the beam are 
zero. This gives four unknown constants in the 
analytical expression for Yrnood - These constants 
are evaluated by applying the four boundary condi- 
tions and the orthogonality properties of the in- 
dividual modes. A normalizing constant called the 
modal mass (Mn) is defined in terms of the mode 
shape and the linear beam mass density P(X). 
Equation (5) shows the relationship between (My), 
(Sn), Wnt) and the driving function £(@). 


y 2 
Mn iy E(x). CO) dx Y 
r 
gr ea {Lec Ga(x) ,0)] dx 
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The exact dimensional values of (Mn) and Y(Xi) 
are unimportant because only the ratio 

LY (xi) YCX;) J Mn] is used in control sys- 
tem transfer functions. When the mode shape is 
determined experimentally, the value of Mn is 
computed in the non-dimensional coordinates de- 
fining n(x). ‘This makes the use of the com- 
puted (Mn) mandatory with the particular modal 
coordinates used. 


Substituting equations (3), (4) and (5) into 
equation (2) by invoking the orthogonality proper- 
ties of the differential equation and taking the 
Laplace transform of equation (2) with respect to 
time, gives the solution for ¥2 (p) which is the 
Laplace transform of Y2(t). ‘The result of these 
analytical steps is given in equation (6) which is 
the transformed total solution. 


(6) 


Go 


Y(x,p)= 2a In Yri(X) & (p) 
N=\ pt + Bip + Wr 


where: F, (p)= Laplace transform of f,(t) 
L 
ree 2% [E00 160) 3 Pn Ca] Pn 
Mn )ax° vOXe 


For a non-uniform bean the integral defining 
(Wn*) is a constant and (Wy) is the scaling 
factor for orthogonality. 


The general solution as given in equation 
(6) for a free-free beam with an experimentally 
derived mode shape is the basis for the applica- 
tion of the Frahm damper. The forcing function 


f(x,t) is derived by solving the differential 
equation of motion of a Frahm damper with a posi- 
tion input and an inertial reaction force on the 
beam as the output. This allows the computation 
of gn and Fo(p) for a particular mode. Due to 
the resonant low pass characteristics of a damped 
spring mass system, the improvement in body modal 
damping can only occur for one mode per Frahm 
damper assembly. Figure 2 shows the system in- 
cluding the input (Z (€)) and output (fe (t) ) 
functions of the Frahm damper. The transfer 
function relating the Laplace transforms of frét) 
to Z (His used in an equivalent feedback system 
to derive the solution for the nth mode. The 
solution provides a means of adjusting the com- 
ponents of the Frahm damper to optimally increase 
the damping ratio of the nth mode. The transfer 
function for fg(t) driven by Z(t) is given in 
equation (7), which also shows the differential 
equations of motion of the Frahm damper of Figure 
ee 


(a) falt) = (Z,-Z)k+BlZ-Z)-mzZ, (7) 


on 22 (Z-2Z,) k+ B(Z- Ze) Differential 
Equations 
(>) Fe(p) = - p*(k+ Bp) 
Zz 2 7 Transfer function 
(P) i +P +k for the reaction 
force 
where: Fp(p) = Laplace transform of fp(t) 
Z(p) = Laplace transform of Z(t) 


The equivalent feedback system of the Frahm 
damper and the missile body is shown in Figure 4 
and only applies to the nth mode. The final 
transfer function relating the lateral beam dis- 
placement to the forcing function for a beam with 
an attached Frahm damper at point (X, ) is given 
in equation (8). 


Sn Fr(p) Yn (x) Lp? + B/m)p + kK/m] (8) 


Nin (x, P) 2 ~ 2 z 

Bip+ Wn |{p*+B p+k J+ [Vou k,8 
[p*+ Biprin JC p +Bp+k] Eevenip [k.B9] 
Application of Root-Locus Theorems 


From Figure (4) the equivalent open-loop 
transfer function GH¢( p) can be derived and is 
given in equation (9). 


(9) 


GH(p)= m Yate *[@/m)p + k/m] 
Mn [p+ B,p+ Wn] [ p*+(2/m)p + k/m] 
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The Root-Locus criteria state that GH(p) 
must equal (-1) at each of the closed-loop roots. 
The variables Y4,(x: ), (Mn), and (B, ) are 
known because the dynamic properties of the mis- 
sile are known and (mM) is the mass of the Frahm 
damper and is known because it is the mass of one 
of the sub-system equipment packages. The values 
of (k ) and (B) must be chosen so that the damp- 
ing of the closed-loop roots have a maximum in- 
crease in magnitude over (B, ) for a given magni- 
tude of (mn). 


The Root-Locus theorems are applied to equa- 
tion 9 along with additional constraints that 
guarantee that the closed-loop roots have a maxi- 
mum increase in damping over the unmodified beam 
equivalent of the missile body. These additional 
constraints are applied to the shape of the Root- 
Locus as (™M) increases from zero and to the lo- 
cation of the open-loop pole and zeros which de- 
pend upon (K/,) and (8/m). Figure 5 shows the 
Root-Locus plot for equation 9. The locus of the 
root leaving the original beam root (Ps, ) is con- 
strained to be an arc of a circle. The maximum 
decrease in (¢) is required for a minimum change 
in (23) to maximize the increase in root damping 
ratio for a fixed (m). This means (dA4¢/dé3 ) 
must be a maximum value. Equation 10 shows this 


evaluation 

Adgd_ | ped Awe dS 

Sil ce Rees Aegean eran, Sa 
where: - 


cos(¢,-4¢) =htor closed-loop roots (R, © R2) 


AR - O for maximum 
Ads (A¢) 


The derivative ( AAp/AL ) has a maxi- 
mum value if (AR/dé3) = 0. Therefore, the 
locus of the root (R, ) must be an arc of a cir- 
cle starting at ( P3; ) with a radius (Wn) as 
shown in Figure 5. To facilitate the analysis, 
the parameters (kK/m) and (B4,) are normalized 
with respect to (Wn). The zero and pole loca- 
tions, and the zero and pole distances, from (R,) 
are given by equation (11) for Figure 5 and for 
equation (10). 


(4) for a given 


Co 
Po f- T+ i veede | o, 
' Vgeee 4 
Lies ee aw@oa te ort 
A> - HaER Jo, 
Pp = [-% +i fa-T*" [Wn 
2 SO Tea 
R= |e i fant” [a, 
iL 2 + 


[ep = (sin (o-aey]"+ [22 — cosg ap} 


] -fes-e re sin(}-a¢)] 4[ 5 -cos6-aff 
ee Aah ae Caen 
] 


S 2 Sila eo) 
[Ea] borer tts sin(}-26)] +[& -c05(¢-ad)} 


where : We =o [kim]? 5 SL res ish - [Bl 


The eight angles shown in Figure 5 are re- 
lated to the variables describing the pole-zero 
lengths of equation (11). Equation (12) shows 
these relationships. 


Ad oa 
= or 
Saree Eb erty 
Tan é€, = sin(¢-ad)—[u,- a a te 
Nef 2 = ees CPHIAD) 
Tan €2 = Tf2 — cos (G-Ad) 
[w.*-@/2) J? + sin(p-A) 
= _Cos(¢-4$) — cos @ _ 
Tan €4 ae 2 ea 
Tanet = s'in(p-sAd) 
W2/T%T —-Cos(¢-49 
fee Bs ee Ades 
pos tee 2 20n beeen = * Zan 


(Damping of closed-loop 
root (R,)) 


ast(Dh — = 
cos(f-Ad) = 5 
The pole-zero lengths and the angles of Fig- 
ure 5 are constrained by the gain and phase cri- 
teria which guarantee that (R,) is a root of 
equation (8). ‘The gain and phase criteria are 


applied to the variables of equations (11), and 
(12) in equation (13). 


(a) 0, = €4+ €5+ €,-€2-2[2-$+44] (13) 


Os Aa 
(b) G= aa ree 


where: (a) Root-Locus phase criterion 


(b) Root-Locus gain criterion 


Ge, ym Pr (xX) 
Mn 


The solution for (7% ), (W.), and ($e) as 
a function of (Go ), (Wn), and (4%, ) is deter- 
mined by substituting the pole-zero lengths of 
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equation (11) and the angles of equation (12) in- 
to the gain and phase criteria of equation (13). 


The solution is not obtainable in closed 
form because the equations for (7 ) and (Wo) 
are not linear, nor do they contain common poly- 
nomial factors which might allow simplification. 


The equations are solvable by a simple con- 
cise iterative cycle. A further optimization on 
the solution for (G% ) and (We) must take place 
during the iterative cycle. The location of (P, ) 
must be approximately the same distance from the 
closed-loop root (R, ) as the original beam root 
(P; ). This means that (49 )~(XL, ). Attaching 
the Frahm damper to the beam introduces an addi- 
tional pair of complex roots (R3) and (R4). 

By constraining the value of (4, ) to be slightly 
larger than (Ls ), and the range of ( €1) to be 
less than (2), the complex roots (R3) and (Rq4) 
will have a slightly larger damping ratio than 
the closed-loop roots (®, ) and (Rz2). These 
eee constraints are exhibited in equation 
14). 


4G Ze Ae (4, approaches £3) 


a (defines the solvable) 
2 (range of the problem) 


(14) 


é, < 


If (€, )> # or (£,)<(4s) then (R3) and 
(Ra) are possibly less damped than (R, ) and 
(R2). This constraint defines the range in (cq) 
and (4d ) in which the optimizing constraints 
hold. The minimum value of ($ -4¢ ) is 
approximately 65° which covers the expected val- 
ues of (4,) and (™). 

Iteration Cycle for (K) and (B 

The iteration cycle begins by choosing a 
value for (4q@) which automatically fixes (£3 ), 
(€,), and (€3) since (¢) and (Wn) are given. 
A value for (€, ) is chosen arbitrarily. ‘The 
value of (€4) is computed from (dp) and (4¢). 
Since (4, ) is approximately equal to £3 , both 
(We) and (Te ) can be determined from the equa- 
tions for (Tan &,) and (4.), which then provides 
the means for computing (€2). The value of (S,) 
is computed from the phase criteria equation (13a) 
The value of ( Wo/ Te )ais computed from the 
known parameters (9, ), sin (6 -4¢) and 
cos (P-4% ). Both (We) and ( Te) are known 
from assuming the value of (€, ). ‘The value of 


(We*/7e )p is computed from (7% ) and (We). A 
comparison is made between (W.*/%s )a and 
(Wo*/% ),. If (Woe /lo )a is greater than 
(wo*/% )% then a second value of (€, ) is 
chosen which is larger than the first value. An 


estimate can be made on the increment in (€, ) by 
using (We*/%e )b to compute a value of (1, ). 
Since (@,,) is larger than (@, ) the difference 
between (6,4 ) and (6, ) is added to the original 


value of (€, ) to start the second iteration 
cycle. The iterative steps per cycle are out- 
lined below. 


Iteration Cycle 


1. Given (¢), fix (Ad) and choose a value of 
(€,) to begin the iteration or take (€ )i 
from the i cycle to commence the (i+!) 
cycle. 

-! [Cos(d-Ad)-Cosd 
€q = Tan Eeorscces| 
As ri sin(4@ ) 
fe 


Wr 


2. Set the ratio Tan(€,)i = [SF] on a slide rule 
and choose specific values of (X%) and (/) 
so that (x *+ B*) ~ 2/wn) 


3. Compute (7. )(: 
Ce)ise 21 Act cos(¢$-A%)] 
4. Compute (Wo he 4 Z 
(jf = [x + sin(d- Ad] +), 
5. Compute (€:2 );: eB 4 
een aes een! 
6. Compute (@,)- 


(Gijinge (64+ (6); += 
7. Compute (Wo /7o Yai? 
s/s = nth con(d- at 
8. Compute ( Wo/ to Yet 
(We/%)b: = (Wo). Atari 
If (We/% )ai>(We/Te ) bi then choose a 


new value of (€, )¢+, larger than (©. )i 
To estimate (€, )i4, : Compute (O,),; 


(6,)4: = Tan"'[-S1n (@- Ad) _] 


(2e2)— Cos(f-4¢) 


eee (€2); “gd 


and 
Cant Fa (E,)¢ + (0,)4.; = (6); 


If (Wo /Zo lai < (Wo / Teo )bi choose a new 
(€,) smaller than (€,){ by step 9, in 
which (@, )4;< ( @)¢ 


10. 


If (We/ )’< {1 then estimate (€, ) c+ 
by setting (€,)-,, = (€); +(@%)p-- (dc 


Za 
With the new value of (€, )i+: re-enter the 
iteration cycle at step 2. 


The completed iteration produces accurate 
values of (T ) and (We), for each pair of the 
fixed values of (}) and (A), which satisfy the 
equations (11), (12) and(13) as well as the opti- 
mizing constraints. The corresponding value of 
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a can then be determined because (“n), CE); 
ki), (4s), (£4 ); bile ), (7), and 
(Y2(x.) Z7Mn ) are known. This relationship is 


given in equation (15). 


Mn 4, £2 Ls £4 
Yh*¢ x.) té L, 


In practice, the value of (m) is given 
along with the other properties of the beam. ‘The 
required inverse functional relationship was 
determined by using the iteration procedure for 
the ranges of (q?) and ( €,) consistent with the 
optimization and solvability constraints and 
plotting the results versus (G,). ‘The result of 
this analytical procedure is given in Figures 6, 
7, and 8. Figure 6 shows the Frahm damper nor- 
malized dashpot coefficient (7 ) as a function 
of the equivalent open-loop gain (Go ) for vari- 
our values of initial beam damping ratio ( %, oe 
Figure 7 shows the Frahm damper normalized spring 
constant (We) as a function of the equivalent 
open-loop gain (Ge) for the range in (4, 
Finally, Figure 8 gives the damping ratio (%, ) 
of the least damped, closed-loop root of the mod- 
ified beam as a function of (Go). 


aan (15) 


Conclusions 


To demonstrate the usefulness of the func- 
tional relationships shown in Figures 6, 7, and 
8, a Frahm damper is applied to a theoretical 
beam whose flexural properties are given in 
Table 1. 


Table 1 
Sample Beam Flexural Properties 
Ya Mode Wn 4, Yn(Nose) Yn( Tail) 
900 slugs 1 30 009 2.0 SS: 
70O slugs 2 [Om O1O -1.0 1.6 


Table 1 could apply to a missile whose lift- 
off weight lies between 75,000 lbs. and 150,000 
lbs. A typical major sub-system module might 
weigh as much as 640 lbs. which corresponds to 
20 slugs. The flexural properties of the re- 
sultant Frahm-damper-modified beam, with the 
damper applied to the nose, is given in Table 2. 


Table 2* 
Flextural Properties of Modified Sample Beam 
(m) (B) (K) 
Mode (4,) (We) (%) slugs lbs sec/ft lbs/ft 
1 .1593 .918 .504 20 303 15,200 
2 .0932 .968 .325 20 456 88, 600 


*(damper applied at nose) 


A comparison of the initial and final values 
of the damping ratio (i.e., (‘,) compared with 
(4, )) indicates a substantial improvement can be 
realized from the Frahm damper principle. The 
spring constant for the first mode as indicated 
in Table 2 is relatively small compared to typi- 
cal lateral spring constants, but would be accept- 
able with its matched dashpot. 


A dashpot with a viscous coefficient of 
500 (lbs. sec/ft) and driven at a frequency of 
4O radians/sec with an input position amplitude 
of 1.00 inches would provide a maximum reaction 
force of 1670 lbs. If the piston area is 1.0 
(sq in), then the maximum internal pressure would 
be 1670 psi with a required piston stroke from (2) 
to (4) inches. A dashpot with these characteris- 


(B)FRAHM DASHPOT 

Uk) FRAHM SPRING 

(m)FRAHM MASS~ \ 
\ 


FRAHM DAMPER 
ASSEMBLY 


(@) Zit 
() f, lt 


+Z, 


! 
i, FRAHM MASS 
-Z, 
@): @ 1S THE COORDINATE MEASURING THE POSITION OF THE FRAHM 


DAMPER BASE 
(b): f- (DIS THE REACTION FORCE EXERTED BY THE FRAHM DAMPER 


BASE DUE TO AN INPUT POSITION (2) at (x=%,) 


Fig. 2. Frahm damper attached to a free-free beam. 
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tics would weigh no more than 10 lbs. Consider- 
ing the potential improvement in beam flexural 
damping this weight penalty would be negligible. 
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Fig. 3. Force function loading of a missile body. 


TRANSFER FUNCTIONS FOR BODY 
FLEXURAL DYNAMICS 


i 
P74 B, P+ QA 


mii) 7°(K+ &e) 
th Bee K) 


Yn (Xs P) 


FRAHM DAMPER TRANSFER FUNCTION 


Fig. 4, Equivalent feedback system for a missile body. 


(+) IMAGINARY 
AXIS 


ROOT-LOCUS OF CLOSED-LOOP 
ROOT EMERGING FROM THE 
ORIGINAL BEAM ROOT(R) 


(~) REAL 
AXIS 


(P,&P,)ARE THE ORIGINAL 
COMPLEX ROOTS OF THE 
UNMODIFIED BEAM 


(R &P,):ARE THE COMPLEX 
POLES INTRODUCED BY THE 
FRAHM DAMPER 


(2,2, €2,)ARE THE REAL AXIS ZEROS 
INTRODUCED BY THE FRAHM DAMPER 3 


(R,& Ro ):ARE THE RESULTANT COMPLEX CLOSED-LOOP 
ROOTS OF THE BEAM 


(Rz& R,)ARE THE COMPLEX CLOSED-LOOP ROOTS 
INTRODUCED BY THE FRAHM DAMPER 


[he] 


Fig. 5. Root-locus of beam with attached Frahm damper. 
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NORMALIZED FRAHM DAMPER VISCOSITY COEFFICIENT (7, ) 


NORMALIZED FRAHM DAMPER SPRING CONSTANT (~,) 


MODIFIED BEAM EQUIVALENT OPEN LOOP GAIN 


Fig. 6. Dashpot coefficient (B) vs open loop gain (Gp) and initial beam damping (C3). 


MODIFIED BEAM EQUIVALENT OPEN LOOP GAIN 


Fig. 7. Frahm spring constant (K) vs open loop gain (G,) and initial beam damping (¢). 
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MODIFIED-BEAM FINAL DAMPING RATIO OF LEAST DAMPED ROOTS 


(4) 


MODIFIED BEAM EQUIVALENT OPEN LOOP GAIN 


Fig. 8. Modified beam final ratio ({,) vs equivalent beam open loop gain (G,) and initial beam damping 
ratio ({). 
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Summa ry 


An analytic study of flywheel autopilots for 
attitude control of space vehicles is described, 
and the results of a three-axis analog computer 
program are presented. 


The objectives of this study were to determine 
the requirements, capabilities and limitations of 
flywheel autopilots as functions of desired ac- 
curacy and speed of response, disturbing mo- 
ments, differential gravity restoring torques, com- 
ponent uncertainties, and vehicle initial attitude 
and attitude rate errors. The stability and ac- 
curacy degradations that resulted from not com- 
pensating for gyroscopic crosscoupling torques 
are summarized. Design tools are presented 
for synthesizing a three-axis autopilot in accord- 
ance with specified design criteria of size, power, 
accuracy, response rate and error disturbances. 


Introduction 


The use of accelerating flywheels for vehicle 
attitude control has been considered and proposed 
for various space missions, particularly when 
long duration flights and accurate vehicle orienta- 
tion are required,.!-3 


The basic equations relating vehicle response 
to flywheel reaction torques and various disturb- 
ing moments have been covered in considerable 
detail. However, the actual synthesizing ofa 
three-axis autopilot system based on various de- 
sign criteria is usually neglected. When flywheel 
controllers are used, the effects of intermode 
coupling between the pitch, rolland yaw axes 
are especially significant. 


This paper presents the results of a study pro- 
gram which was undertaken to provide a better 
understanding of the system synthesis problem, 
including the effects of intermode coupling. 


Control Application 


The ability of inertia wheels to control vehicle 
attitude stems from the basic law of action and 


reaction. If an inertia wheel mounted within the 
vehicle is accelerated about its axis of rotation, 
the vehicle is accelerated in inertial space about 
the same axis, but in the opposite direction. 
Assuming zero initial conditions, 
‘6 = -Io/P (refer to list of symbols). (1) 
Utilization of this basic principle to simultane- 
ously control vehicle attitude about all three axes 
requires an autopilot consisting of attitude dis- 
placement sensors and/or rate sensors, control 
wheels, and electronic circuitry which must 
govern the rotation of the wheels as a function of 
the attitude displacement and/or rate signals. 
Figure 1 presents a functional block diagram of a 
single-axis control loop. Three such loops are 
ordinarily required, one for each axis. The at- 
titude displacement sensor could be a gyroscope 
or an inertial platform for short duration missions, 
but in long duration missions (for which flywheels 
are especially attractive) gyro drift rate would 
have to be corrected by means of some monitor - 
ing device such as an horizon scanner or a stellar 
tracker. Aerodynamic forces are negligible when 
traveling in space and, consequently, attitude 
damping must be artificially provided by use of 
rate signals, either from a rate gyro as shown in 
Fig. 1 or from a lead circuit operating on the 
displacement error signal. Flywheel accelera- 
tion, in response to the sensor error signals, pro- 
vides the restoring torque necessary to control ve- 
hicle attitude in the presence of initial condition 
errors and disturbing moments. Gyroscopic 
torques about each axis result from rotations of 
the vehicle and of internal machinery about the 
other two axes. 


From the viewpoint of system weight and 
complexity, the differential gravity torque shown 
in Fig. 1 provides an attractive method for sta- 
bilizing an earth satellite to the local vertical (in 
a low altitude near-circular orbit). This method 
will be subsequently discussed in some detail. 


The coordinate reference frame employed and 
the desired orientation of the vehicle are impor- 
tant considerations, Figure 1 is applicable with 
either an earth-fixed local vertical system or an 
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inertial frame fixed with respect to the stars. 
The former is desirable for a communication or 
reconnaissance satellite, whereas the latter is 
more advantageous for an astronomical satellite 
or an interplanetary space probe. If a local 
vertical orientation is maintained, additional 
crosscoupling terms are introduced as a result 
of the rotation of the coordinate reference frame 
in inertial space. 


Study Objectives and Approach 


In order to evaluate the requirements, capabil- 
ities and limitations of flywheel controllers, an 
analytic program with the following objectives 
was undertaken, 


(1) What is the most efficient form of the 
controller transfer function (ratio of flywheel 
angular velocity to sensed error signals)? 


(2) What compromises are possible as 
functions of attitude accuracy requirements, 
speed of response, initial condition and compo- 
nent inaccuracies, disturbing moments, wheel 
saturation, and system size, weight and required 
power? 


(3) What stability and accuracy degrada- 
tions result from the gyroscopic crosscoupling 
torques ? 


(4) What is the best way to compensate for 
these crosscoupling moments and what complexity 
is involved? 


The analytic study was accomplished by means 
of a two-phase program. Phase I involved single- 
axis hand analyses for approximately determin- 
ing the first two objectives, whereas a three-axis 
analog computer program was undertaken to 
evaluate the crosscoupling effects. 


In view of the study objectives noted, no 
attempt was made to determine an optimized atti- 
tude displacement error sensor. Depending upon 
the specific space mission, various types would 
be applicable and it was presumed that any neces- 
sary coordinate conversions could be accomplish- 
ed with zero error, 
on two types of attitude reference systems. The 
first assumed a desired orientation fixed in 
inertial space and a continuous source of three- 
axis, attitude displacement error signals with- 
in the space vehicle. The second assumed a 
coordinate reference frame defined by the earth's 
local vertical and orbital plane (that is, rotating 
in inertial space) but with no source of continuous 
displacement signals available in the vehicle, 

The latter would utilize differential gravity as a 
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The investigation concentrated 


restoring torque, with rate gyros and flywheels 
for attitude damping. For this orientation, the 
programmed orbital angular velocity would have 
to be subtracted from the rate gyro output signal. 


Single-Axis Investigation 


The optimized form of the controller transfer 
function was determined by studying the transient 
response of the system to various forcing func- 
tions and initial conditions. Referring to Fig. 1 
and assuming a displacement error sensor, the 
idealized controller transfer function is 


hed ete (2) 


The integration 1/S is needed to prevent the 
attitude error (9) from building up due toa 


small but constant disturbing torque (M,)- For 
example, if o/e , has the form -1/I instead of 


-1/1S, then 


g -1 


Me tS Se Sie ye 
= 8 fi“ | 


: 1 


(3) 


Another advantage of using the controller de- 
fined by Equation (2) is its suitability to cross- 
coupling torque compensation. For example, 
if the gyroscopic moments are measured and 
inserted as compensating signals into the con- 
troller (Fig. 1), the resulting control moment 
(Md) would exactly equalize the disturbing 


torque. 


The characteristic equation for the one-axis 
attitude control loop, assuming Equation (2) and 
a displacement error sensor, is 

us 2 


ares S+ £- s°=0. (4) 
1 1 


K 


Speed of response and accuracy are both im- 
proved as K, is increased and system damping 


is governed by K. For example, a damping 
ratio of 0.5 is attained by making 


1 
Soe (K, P) oh (5) 


Transfer functions of system response to vari- 
ous forcing functions are shown in Table 1, and 
transient response plots (for example, Fig. 2) 
for each equation in Table 1 were determined 
for various values of the significant parameters. 
Table 1 indicates that the steady-state value of 


0. that would result from a constant disturbing 


moment or rate gyro drift would vary inversely 
with displacement gain (K,). The primary 


problem associated with this control scheme is 
that the flywheel continuously accelerates in the 
presence of a constant disturbing moment (as 
noted from Equation (4) of Table 1). There 


appears to be no controller form which can remedy 


this situation, unless some additional restoring 
torque (such as differential gravity, solar sail 
or aerodynamic force) is introduced. Conse- 


quently, these disturbing torques, although small, 


would have to be equalized for this particular 
type of system by fixed reaction jet impulses 


which would be fired when the flywheel approaches 


saturation, 


The use of differential gravity to remedy fly- 
wheel saturation was investigated. For this 
type of system which utilizes only a rate gyro 
and flywheel controller, the optimum form of 
the controller was determined to be 


-K 
ey Vee eae ; 
“Ko 
For this case, a pure integration such as are 


although more desirable from the viewpoint of 
system simplicity and stability, would cause the 
flywheel to continuously accelerate due to a rate 
gyro drift or an error in the preprogrammed 
orbital rate. 


The characteristic equation, using Equation 
(5) with Ko = 7, /I, is 


(P + - K.) 
Ce Ss + s*+ CBE 2 2G! 
4 4 


where the differential gravity restoring torque, 
as derived for the small angular displacements 
from the local vertical*, is defined for the 
pitch (y) axis as 

(7) 


2 
eo 
Redrawing Fig. 1 as Fig. 3 simplifies the sys-~- 
tem synthesis problem, which is essentially the 
selection of optimum values of ; and K,- P and 


K {P? = PF). 
x Zz 


K, are fixed by the vehicle configuration and 


orbital period (a near-circular orbit is required 


for this type of system to be effective). Referring 


to Fig. 3, satisfactory system damping is pro- 
vided by making 
Perks 
2~ 2 = 
Tv = = ae and T Ky = 


1p, 


249 


Transfer functions defining the response of 
such a system to various forcing functions are 
shown in Table 2,and Fig. 4 shows a typical tran- 
sient response plot of vehicle attitude error and 
flywheel momentum. Table 2 indicates that the 
flywheel will not continuously accelerate in the 
presence of a steady-state disturbing torque 
(Fig. 4) or rate gyro drift. However, because 
of the small magnitude of Ky as defined by 


Equation (7), system accuracy and speed of re- 
sponse are, of necessity, quite low. A conse- 
quent limitation to the use of this type system is 
that initial conditions (in particular vehicle 
angular rate) must be fairly precise. As will be 
subsequently described, this problem is in- 
creased by an order of magnitude for the practi- 
cal case of simultaneous three-axis control. 


Three-Dimensional Equations of Motion 


The three-dimensional equations of motion 
for the case in which the desired attitude 
orientation remains fixed in inertial space, as 
derived in Appendix A, are presented in Table 3. 
It is noteworthy that the preceding single -axis 
analyses neglected the final three terms on the 
left side of each equation. These terms repre- 
sent the gyroscopic crosscoupling torques--the 
product of angular momentum about one axis 
multiplied by angular velocity about the second 
axis produces a torque about the third axis. 


For the desired orientation defined by the 
earth's local vertical and orbital plane, the sys- 
tem equations become considerably more com- 
plex. However, by assuming small angular 
deviations from the desired vehicle orientation 
and a constant orbital angular rate (wo), the 


equations can be managed. These are derived 
by making the following substitutions into the 
equation of Table 3. 


8 = 6) - % B= 6, 
Sa Pte Vu © $= 4) > wy Hy, 


A further simplification (a, = 0) is attained by 


using a minimum equipment system? employing 
two flywheels and two rate gyros. The rate 
gyros sense angular velocities about the vehicle 
pitch and roll principal inertia axes, whereas 
the flywheels are aligned along the pitch and yaw 
axes--the pitch axis is defined as being normal 
to the orbital plane and the yaw axis through the 
center of the earth for the case of zero attitude 
error. With this configuration, roll control re- 


sults from gyroscopic coupling (through wo) 


with the yaw control mode which does use a fly- 
wheel. The differential gravity restoring torques, 
which must be added to the equations of Table 3, 
are as follows. 


2 
M, = 309 (Py Ee core 

x 
M essa -P)e 
E e0- x ee & 

y 

=0 9 
Mp (9) 


Three-Dimensional Simulation Program 


The gyroscopic crosscoupling effects were 
investigated by mechanizing the three equations 
of motion in The Martin Company analog-com- 
puter facility. This program utilized the con- 
troller transfer functions which were optimized 
by means of the single-axis synthesis, with 
allowances for gain and time constant variations. 
Provisions were made for varying all of the 
parameters shown in Fig. 1 and for introducing 
various inputs into each of the three control 
loops, either individually or simultaneously. 


The gyroscopic moments, computed in accord- 
ance with the equations in Table 3, were intro- 
duced as disturbing moments and also, in some 
cases, as compensating signals to the flywheel 
controllers. In addition to constant disturbing 
torques, impulse torques were introduced to 
simulate a meteoritic impact, an initial velocity 
error or a change in angular momentum asso- 
ciated with vehicle rotating equipment. 


A summation of the significant results of 
this simulation program follows. In addition to 
the analog study, the roots of the roll-yaw 
characteristic equations were determined by 
digital computation for the system employing 
differential gravity and two wheels. These re- 
sults showed excellent correlation with the 
response plots from the three-axis analog pro- 
gram. 


Computer Results of the Three-Wheel System 


The three-axis analog simulation of the sys- 
tem utilizing an inertial frame of reference and 
three flywheels fixed to the vehicle principal axes 
revealed little or no stability degradation due to 
the crosscoupling terms, with or without com- 
pensation. However, the gyroscopic moments 
(when not compensated) caused significant tran- 
sient errors in vehicle attitude, but in all cases 
these errors subsequently damped to zero. Sys- 
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tem response to attitude displacement com- 
mands, disturbing torques and rate gyro errors 
was determined for a wide variation of gain set- 
tings--system period of oscillation ranged from 
60 to 2800 seconds for all axes. The transient 
response plots were as expected from the single 
axis analyses (Table 1). 


A typical transient response plot is shown in 
Fig. 5, for which the following set of parameters 
was used. 


= 1000 slug-ft” 


Pa 
x ¥ 
2 
P. = 100 slug-ft 
(@ t inh ( h ‘ ): Os 4 ee ee 
ontroller (each axis): = = 3a +0.58) 
K, was fixed so that the approximate natural 


frequency was 0.1 radian/second for each axis. 


Ky was fixed so that the approximate damping 


ratio was 0.5 for all axes. 


Figure 5 plots system response to two types 
of inputs. The first, at time t = 100 seconds, 
was a pitch impulse disturbance of 11 foot- 
pound-seconds (equivalent to a vehicle initial 
angular rate of 0.63 degree/second in inertial 
space) and assumed net angular momentums of 
zero about the roll and yaw axes. As expected 
from the single-axis studies (Table 1), the pitch 
flywheel absorbs this momentum change after 
which pitch attitude is stabilized to zero error. 
From this and similar transient plots, a good 
approximation for determining peak overshoot, 
assuming a 0.5 damping ratio, is 


-AH 
2(PK yi 
1 


(0 ) = 


b (10) 
peak 


/2 


where AH = M or AH = PO. If the damping 


D/S 
ratio is decreased by using a lower value of K 


( 0 ) would be increased. 
peak 


Ae 


System response to the second input, a roll 
displacement command of 10 degrees, illustrates 
the gyroscopic crosscoupling effects. The peak 
yaw displacement of 0.084 radian in Fig. 5 re- 
sults from the induced yawing moment defined 
by - Le From Fig. 2. ~ - 11 foot pound 
seconds and the average ¢ for a 30-second time 
interval is approximately 0.007 radian/second.. 
Substituting these values into Equation (10), 


(p ) is determined to be 0.082 radian, which 


« peak 


explains the observed results. Gyroscopic 
coupling into the pitch axis is seen to be consider - 
ably less because C ae Re (Equation (10)) is 
greater by an order of magnitude than Wee Rg, 
The residual yaw wheel momentum of 2 ft-lb-sec 
results from the gyroscopic torque (¢ ve ) 


applied for 30 seconds in one direction. This 
increment of yaw momentum was absorbed by the 
yaw flywheel. When roll displacement was com- 
manded back to zero, yaw wheel momentum re- 
turned to zero. 


When gyroscopic torque compensation was 
applied, the response of all three flywheels, as 
plotted in Fig. 5, remained essentially the same. 
However, the yaw attitude excursions were elim- 
inated by introducing these compensating sig- 
nals (the final three terms on the lett side of 
each equation in Table 3) into the flywheel con- 
trollers. In other analog runs with the controller 
gain set atl. 25/1 instead of 1.0/I, compensa- 
tion was approximately 75% effective. 


Simultaneous attitude displacement commands 
of 10 degrees about each axis were also intro- 
duced. Response was quite satisfactory even 
without gyroscopic torque compensation; some 
degradation in transient response resulted, 
particularly when an initial angular momentum 
was present, but all axes stabilized to zero with- 
in the expected time duration. A typical re- 
sponse for a low gain system is reproduced in 
Fig. 6. The noticeable lag in yaw response is 
attributed to the comparatively low vehicle 
inertia about the yaw axis together with the gyro- 
scopic torque ( d ie a). With this inertial frame 


of reference, the recommended vehicle configura- 
tion would have P. = i = Pi This inertia 


ratio is almost mandatory in a low-orbit earth 
satellite because of differential gravity torques. 
The alternative would be to utilize differential 
gravity to help provide a local vertical orienta- 


tion with os = Py => Pi. 


Computer Results of the Two-Wheel System 


Initial analog results of the two-wheel sys- 
tem, using differential gravity together with 
pitch and yaw controller transfer functions as de- 
fined by Equation (5), indicated marginal stabil- 
ity. A roll-yaw damping ratio of approximately 
0.15 was the best that could be obtained by vary- 
ing + and Ky, whereas pitch response was as ex- 


pected from the single-axis analyses (Fig. 4). 
Because of this roll-yaw stability degradation 
and the increased complexity introduced into the 
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analog simulation by Equations (8) and (9), the 
analog program was terminated after the initial 
runs to conduct additional analysis. The three- 
axis system, as simulated, was reduced toa 
two-axis roll-yaw coupled mode by assuming o 


and @ constant. The roots of the resulting sixth 
order characteristic equation were extracted 
and, in all cases, good correlation (in terms of 
frequency of oscillation and speed of convergence 
or divergence) was established with the analog 
response plots. The correlation demonstrated 
the effectiveness of this root extraction method 
for system optimization and also substantiated 
that a 0.15 damping ratio was all that could be 
expected. 


The form of the roll-yaw characteristic 
equation is developed as Equation (B-3) in 
Appendix B. Tables 4 and 5 illustrate the re- 
sulting roots for some typical cases and also 
show the criticality of ; and K,. Varying these 


values by a factor of two can change an 
optimum damped system to one of negligible 
damping or actual instability. 


Table 5 shows that stability need not be de- 
creased by incorporating a constant wheel mo- 
mentum on the pitch axis, although fluctuations 
of the yaw control parameters 7, and Ky (which 


must be selected for this constant value of 7 es 


appear to be even more critical. This increased 
pitch momentum, having the same polarity as 
Wys provides additional yaw stiffness and thereby 


reduces the yaw transient error resulting from 
an initial yaw or roll rate or a subsequent dis- 
turbing torque. The analog response plots veri- 
fied that the peak overshoot in yaw could be 
approximated from Equation (10). The actual 
peak values were somewhat greater due to damp- 
ing of less than 0.5. In applying this equation, 


K, must be replaced by the differential gravity 
restoring torque (K,) about the yaw axis, which 


can be approximated (see Appendix B) as 

ws Vervieea 
where Wp ve is equivalent to the constant value of 
pitch wheel momentum Bs os Thus, witha 
constant i eee of -0.9 ft-lb-sec as assumed in 


Table 5, it follows from Equation (10) that the 
peak yaw overshoot would be decreased by 2. 24 
as compared to the case with zero vA 


Physical Characteristics 


The physical requirements of the system in 


terms of flywheel size, weight and power depend 
upon the maximum wheel momentum needed to- 
gether with its rate of change. From Equation 
(1), the wheel momentum must be sufficient to 
provide a specified maximum slewing rate of the 
vehicle. However, if vehicle stabilization alone 
is required, the wheel need only equalize dis- 
turbing torques and initial conditions. For this 
situation, the flywheel physical requirements can 
be held to a minimum by compromising system 
speed of response, K, should be just sufficient 


to provide the specified static accuracy. The 
equations in Tables 1 and 2 together with accurate 
predictions of initial conditions and disturbing 
torques provide the basis for selecting minimum 
size components consistent with the vehicle con- 
figuration and desired performance. 


A parametric plot of flywheel size and weight 
is presented in Fig. 7 as a function of required 
angular momentum. Figure 8 shows required 
power as a function of wheel momentum and its 
rate of change where 


ile HE ie (NG)- 


Watts (average) = 5t 


eegiseall ‘Cssouie 


10°t 


(12) 


and t is the time required to accomplish a Ao 
(radians/second) change in wheel velocity. 
Minimum power is obtained by using a maximum 
inertial wheel, as indicated from Equation (12). 


Typical values of response time and flywheel 
momentum are shown in Figs. 2, 4, 5 and 6. 
These values illustrate the significance of sys- 
tem gain and vehicle inertia in determining re- 
quired physical characteristics. 


The recommended controller transfer func - 
tion defining the ratio of flywheel response to 
sensed error signals (a/« BEA 1/IS) can be 


realized using off-the-shelf components. A 
direct-drive motor and electronic circuitry in- 
corporating a single integration would suffice for 
each axis. Circuit design and repeatability are 
not critical. A controller smoothing time of as 
much as 0.5 second and a 25% gain variation 
were shown to be tolerable. 


If dynamic accuracy requirements necessitate 
compensation for gyroscopic moments, rate 
generators fixed to each flywheel shaft and the 
associated computing circuitry would have to be 
incorporated. 


Applications and Problem Areas 


The single-axis analyses and computer studies 
previously described made use of certain lineariz- 
ing assumptions (such as small angle approxima - 
tions). Also, various second order disturbances 
(such as orbit eccentricity) were neglected. 
Eccentricity effects could be serious with a low 
gain system and, in particular, when differential 
gravity is used; the control system natural fre- 
quency would be of the same order of magnitude 
as that of the sinusoidal input disturbance. 


Although more complex computer programs 
are certainly recommended, it is believed that 
such effects can reasonably be inferred from the 


’ 


results of this and other publications. For 
example, computer transient response plots 


of vehicle oscillation about the local eral 
for an undamped "dumbbell" configuration in- 
dicate half amplitude oscillations of 6 degrees 
peak for a 7% eccentric orbit. Moreover, it 
appears that the amplitude of this oscillation 
would not be substantially reduced by the fly- 
wheel damping system--computer response plots 


: : 6 ; 
of a mechanical passive damper substantiate 
this conclusion. 


As an outgrowth of the analytic results pre- 
sented, it seems pertinent to point out various 
space missions for which each of the two types 
of flywheel controllers are recommended. If 
precise vehicle orientation is required, such as 
in astronomical satellites or during periods of 
stellar navigation ''fixes, " the controller form 
defined by Equation (2) should be employed. This 
permits the use of a relatively high gain system 
to minimize errors in accordance with the equa- 
tions of Table 1. It also provides a low-torque 
linear control capability. However, for each 
particular application requiring precise accuracy, 
the actual controller hardware should be bread- 
boarded and included in a three -axis analog com- 
puter study to evaluate the effects of any existing 
nonlinearities. The predominant drawback is the 
tendency of the flywheels to saturate. This prob- 
lem could be remedied by using fixed low-thrust 
reaction jets to desaturate the wheels (such jets 
would probably be present for initial angular rate 
stabilization),If this controller is employed to pro- 
vide a local vertical orientation, such as witha 
high-accuracy reconnaissance satellite, the re- 
lationships in Equations (8) and (9) would have to 
be included in the three-axis equations of motion 
as presented in Table 3. Although this degree 
of sophistication was not achieved in the analog 
program, no stability or accuracy degradations 
are to be expected unless an extremely low dis- 
placement gain is used. 
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If attitude deviations from the local vertical 
of between 1 and 10 degrees are permissible, as 
might be the case with a communication or 
meteorological satellite, the minimum equip- 
ment system utilizing differential gravity would 
have application. The recommended form of the 
pitch and yaw controllers is as defined in 
Equation (5). As noted from Table 2, the satura- 
tion problem is circumvented, but other prob- 
lems are introduced because of the very low 
static gain (differential gravity restoring torque). 
Consequently, the effective use of this type con- 
trol presupposes extremely low initial angular 
rates and a near-circular orbit. Another attrac- 
tive use of this type system, as in a reconnais- 
sance satellite, would be as a backup in the 
event the airborne attitude reference system 
should malfunction. In addition to simplicity of 
operation, the desired vehicle configuration 
es Fo >> P) provides the possibility of 


indefinitely maintaining a local vertical orienta- 
tion after all electronic equipment has been shut 
down. One other "dumbbell" application would 
be for crude alignment of an earth satellite, with 
fine alignment provided by gimbaling the recon- 
naissance sensors. 


Conclusions 


Principal results and accomplishments are 
summarized as follows. 


(1) Three-axis flywheel control of vehicle 
attitude can be effectively employed. The low- 
torque linear control characteristics provide an 
excellent accuracy capability. 


(2) When the desired vehicle orientation is 
fixed in inertial space, gyroscopic coupling has 
negligible effect on system stability. Compensa- 
tion for such effects can significantly improve 
dynamic accuracy, but the additional complexity 
is not warranted for most missions. 


(3) With a desired local vertical orienta- 
tion and a low response control system (that is, 
differential gravity), stability is significantly 
affected by gyroscopic coupling. 


(4) Effective design tools have been devel- 
oped, with illustrative examples, for synthesiz - 
ing a three-axis flywheel autopilot as functions 
of desired accuracy and speed of response, 
vehicle inertia characteristics, disturbing torques 
component inaccuracies, vehicle initial attitude 
and attitude rate errors, and flywheel size, 
weight and power requirements. 


(5) Major problem areas are system reli- 
ability during long missions and flywheel satura- 
tion. 
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(6) With a low-altitude near-circular orbit 
the use of differential gravity minimizes both of 
these problems, provided precise initial condi- 
tions can be achieved. 
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List of Symbols 


I Flywheel moment of inertia 

o Flywheel angular velocity 

ip Moment of inertia about principal 
inertia axis of vehicle 

4 Principal inertia axes of vehicle, 


assumed orthagonal in this study 
o Vehicle rotation (in inertial space) 
about x axis 
Vehicle rotation (in inertial space) 
about y axis 
Vehicle rotation (in inertial space) 
about z axis 
Flywheel reaction torque acting upon 
vehicle 
External torque acting upon vehicle 


fe) 


Q 


ies) 


Total moment acting upon vehicle 


4g 


Angular momentum 
Laplace operator 
Gain factor 

Time constant 

Time or time interval 
Vector cross product 


Seana eS SS 


Wy Orbital angular velocity 

Wy, Vehicle angular velocity in inertial space 

if Equivalent value of constant wheel inertia 

y along the y axis 

4 Error signal to flywheel controller 

a Attenuation factor 

A Increment of 

rpm Revolutions per minute 

G Rate gyro inaccuracy or error in prepro- 
grammed angular velocity of vehicle 

C(S) System characteristic equation 

R Flywheel radius of gyration 

Subscripts 

L Refers to vehicle motion with respect to 
local vertical orbital-plane coordinate 
system defined (for zero error) as hav- 
ing the z axis through the earth's center 
and the y axis normal to the orbital 
plane. 

F Refers to flywheels 

V Refers to vehicle 

c Refers to command (or zero reference) 
attitude 

€ Refers to attitude error 

il Refers to initial condition 

u Refers to unit vector 

x, y,zZ Refers to principal inertia axes of 
vehicle 

Superscripts 


Vector quantity 
d/dt 
d7/at” 


TABLE 1 


One-Axis Transfer Functions with Inertial Co- 


ordinate System 


1 
Cc (S) ky) 
6 
€ -1 
a, PS K, C(S) (2) 
K 
2} 
K,¢c® (3) 
eae sisi t ghee = (4) 
6. PS S C(S) 
£ 
K, PS 
K, cis) 8) 
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Controller Form: 


(1) 


(2) 


(3) 


(4) 


cole ee i 
€ 1 mi Is 
K 
C(S) see BEDE, os 
1 1 
TABLE 2 
Pitch-Axis Transfer Functions with Differential 
Gravity 
oii hey wich 2 wa mL Annas) 
M oA. PS K, C (S) 
i 
a : K, TS 
G K, C(S) 
io: Ia: K S 
Ve ean 20 Ries 
Ma Or Ps K, C(S) 
1 
Io -¢ K, (1 +P _/K s*) 
ide 2 5 
G C(S) 


Controller Form: 


1/2 
a le Se hens ay 
af debe ae K, 
ECG -P) 
4 einen Zz 
(qK, +P) Toke 
OO wn RAIDER Ds AOSTA 
ei Lo 


TABLE 3 


Three-Dimensional Equations of Motion 


I e a P ee i e x . e e _ 2 A 
x o. = d + 90 oO. 1s w he & +0 (P. a) My (roll axis) 
x 


lvoe oP Oe oe ad +h - : 
ae - Me e oo, Ls pias Es, P) M (pitch axis) 


[Nor eb marc hie Ola.1l. 400 6 Beery ie 
7 oy : uD ae es Gale b ee P) M,, (yaw axis) 


I = Flywheel inertia 

o = Flywheel angular velocity 

8,6, = Vehicle pitch, roll and yaw inertial rates 
P = Vehicle inertia 

Mp = External torques 


x, y and z refer to the roll, pitch and yaw principal inertia axes. 


TABLE 4 


Roots of Roll-Yaw Characteristic Equation with Me =0 


K Minimum 
1 (sec) aes Roots Damping Ratio 
100 2225 -0.0081,  -0.010, -0.00017+ 0.0021, -0.0007 + 0,0010 0.081 
200 oan -0.00055, -0.0050,  -0.00010+ 0.0021, -0.0021 + 0.0028 0,048 
400 2.25 -0.00023, -0.0028,  +0.000062+ 0.0019, -0.0011 + 0.0039 Divergent 
-0,033 
400 4.5 -0.00011, -0.0025,  -0.000012+0.0021, -0,0012 + 0.0053 0.006 
400 L125 -0.00047, -0.0025, -0,00015+0.0021, -0,00085 + 0.0025 0.071 
200 1,125 -0.0031, -0.0050, -0.000094+9,0022, -0.00086 + 0.0011 0.043 
400 0.5625 -0.0011, -0.0025, -0.00015+90,0022, -0.00055 + 0.0015 0,068 
200 0.5625  -0.0042, -0.0051,  -0.000076+ 90.0022, -0,.00031 + 0.0012 0,034 
50 0.5626 -0.0189 + 0.0045,  -0.000019+ 0.0022, -0,00007 + 9.0011 0.009 


2 
< 2000 slug ft, Py = 1800 slug-ft-, Ee = 200 slug-ft 


w. = 0.0011 rad/sec 


0 
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Roots of Roll-Yaw Characteristic Equation with at Wo = 


TABLE 5 


-0.9 Himeic ese 


Minimum Damping 


zt (sec) Bo Roots Ratio _ 
100 Hh dhs) -0.0092 + 0.0029, -9. 000083 + 9.0031, -9.00020 + 0,0018 C2027 
200 1.125 -0.0042 + 0.0017, -9.00050 + 9.0028, -0.00033 + 0.0022 0.15 
400 25 -9.0017 + 0.00062, +0,0000004 + 0.0020, -0.00085 + 9.0040 Divergent 
-0.0002 
400 0.5625 -0.00041 + 0.0035, -0,000040 + 0.0019, -0,0021 + 0.0007 0,021 
200 0.5625 -9.0049 + 0.0017, +0, 000033 + 0.0034, -9.00016 + 0.0016 Divergent 
-0.010 
200 25 20 -0.0021, -0.0050, -0. 000028 + 0.0018, -0.0014 + 0, 0046 0.016 
800 2.25 -0.00042, -0.00125, -0, 000054 + 0.0019, -0,00042 + 0.0054 0.028 
800 4.5 =07 00025, -05 0019, -9.00091 + 0.0021, +0,00070 + 0,0046 Divergent 
-0.015 
400 4.5 -0.00051, -0.0025, -9.000026 + 0.0018, -0,00098 + 0.0068 0.014 
200 4.5 -0.0012 + 0.00045, +0,00037 + 0.0016, -0.0042 + 0.0126 Divergent 
=0, 23 
2 2 2 
a = 2000 slug-ft , ee = 1800 slug-ft , Pe = 200 slug-ft 
W) = 0.0011 rad/sec 
APPENDIX A eX = = y., = Oz. (A -6) 
THREE-DIMENSIONAL EQUATIONS OF MOTION wy, x Nea oz - Shy (A-7) 
The moment equations for a space vehicle o,, X Zz = Ox, = oy, (A -8) 


using three variable-speed flywheels fixed to 

the vehicle principal inertia axes, assuming con- 
stant vehicle and flywheel moments of inertia, 
are derived from the following vector equations 
(refer to list of symbols). 


My = H = Hi, + Hy (assumes M,, = 9) 

a = Ilo 

ce =s Oy 

eee yo, z,,*0, (6, xx.) 
toy (a Xy,) aor (0 Xz) 


Vv $x, + oye Pee + (0 Xx ) 


+0 (oy Xy,) + (OX z) 


(A-1) 
(A-2) 


(A-3) 


(A-4) 


(A-5) 
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Equation (A-5) assumes a coordinate axis sys- 
tem fixed in inertial space. 


Combining the above equations according to 
the x, y and z components gives the complete 
moment equations about each of the three princi- 
pal inertia axes. These equations are presented 
in Table 3. The external torques (M,,) would 


arise from such sources as atmosphere, solar 
radiation, mass attraction, etc. 


APPENDIX B 


ROLL-YAW CHARACTERISTIC EQUATION WITH 
ROTATING COORDINATE AXIS SYSTEM AND 
DIFF FERENTIAL GRAVITY 


The equations shown in Table 3 must be modi- 
fied if the desired orientation of the vehicle is 
rotating in inertial space. An important applica- 
tion is an earth satellite, the body axes of which 
are to remain aligned with the local vertical. If 

9 is assumed constant, the necessary modifica- 


tions are defined by Equation (8). Also, differ- 
ential gravity torques must be considered unless 
‘y Po For small displacements from 


160) 


the desired orientation (z axis through the 
earth's center and y axis normal to the orbital 
plane), the gravity torques are defined by Equa- 
tion (9). 


For a ''minimum equipment" system requir- 
ing only two control wheels and two rate gyros, 
the roll-yaw coupling (through %) can be effec- 


tively analyzed by assuming constant and 


Je 
zero @, 


Using these approximations and Equations (8) 
and (9) together with 


o ~ 0 
56 
thio te iw 
yy 0 
he mala 7 Ky (7 = «i 4) 
°, fe Ces) 4 


the first and third equations in Table 3 reduce to 


ek = =~ cpap 
y, |. (P Bs ee ae 
29 Ky t 
Sg arsS iy 
@ K T 
2 OL 2 2 
—- St =e er) |= 
Yr [P,s a SR 2 “0 ee x »| 
K, x 
2 2 
+ P -P -P -I')S 
lees § Oe can ty: ode) 


(B-2) 
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By substituting (B-2) into (B-1) and expanding, 
the roll-yaw characteristic equation is determined 
to be 


6 2 
SEM + S.2P" Pp oe 
SFA Zi 


+s*[P. at wyt (9 {4P,K,+P_K +KoI 
S49) feeape Kt) 


3 2 
+ 
S* wyt [eo @P, Ket 8P K+ 2K, ) 


6 


+ - 
Ky (Fo a P.. K)| 


+4y,K 


oer 


K we ne, : 


2S wo [oot (Ky | 6 


2 
+ + + 
ce SS a us 5 ] 


3 
+ - + - 
Swot [ Bay Kg Ks K, (K, aK, K,)| 


+400 K <. 20 (B-3) 


5) 


where 


This sixth order equation can be used as a de- 
sign tool for optimizing Ky, qt and i These 


are the critical parameters since the others are 
fixed by the chosen orbit and vehicle configura - 
tion. Again it must be emphasized that the orbit 
should be near circular and that ae P>> ss 


for optimum performance. A straightforward 
method of determining system stability involves 
the substitution of parametric values into Equa- 
tion (B-3) and subsequent root extraction or 
application of Routh's stability criterion. Al- 
though this is a tedious procedure, it was ac- 
complished for selected values of Ky, 7 and Ly 


in order to evaluate their significance. Tables 4 
and 5 present the roots, as extracted by digital 
computation, for various combinations of these 
critical parameters. These roots demonstrate 
that the synthesis tools developed during the 
single-axis hand analyses are helpful for "'ball- 
park" selection of parameter values, but that sys- 


tem stability should be verified by the use of 
Equation (B-3). 
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